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NONUNIFORM  SHEAR  STRAINS  IN  TORSIONAL  KOLSKY  BAR  TESTS 

ON  SOFT  SPECIMENS 

Adam  Sokolow  and  Mike  Scheidler 

We  investigate  inertial  effects  in  torsional  Kolsky  bar  tests  on  nearly  incompressible,  soft  materials.  The 
results  are  relevant  for  materials  with  instantaneous  elastic  shear  modulus  on  the  order  of  1-1000  kPa 
and  density  on  the  order  of  water.  Examples  include  brain  tissue  and  many  other  soft  tissues  and  tissue 
surrogates.  We  have  conducted  one-  and  three-dimensional  analyses  and  simulations  to  understand  the 
stress  and  strain  states  that  exist  in  these  materials  in  a  torsional  Kolsky  bar  test.  We  demonstrate  that 
the  short  loading  pulses  typically  used  for  high  strain-rate  (e.g.,  700/s)  tests  do  not  allow  the  softer 
specimens  to  “ring-up”  to  uniform  stress  and  strain  states  and  that  consequently  the  shear  stress  versus 
shear  strain  data  reported  in  the  literature  are  erroneous.  We  also  show  that  normal  stress  components, 
which  are  present  even  in  quasistatic  torsion  of  nonlinear  elastic  materials,  can  be  amplified  in  dynamic 
torsion  tests  on  soft  materials. 


1.  Introduction 

Kolsky  bar  tests  are  widely  used  to  study  strain-rate  effects  in  inelastic  solids  [Gray  2000;  Chen  and  Song 
2011].  In  a  compression  Kolsky  bar  test,  a  relatively  thin  specimen  is  sandwiched  between  two  bars. 
Impact  at  one  end  of  the  incident  bar  generates  a  compressive  wave  that  travels  along  the  bar.  through 
the  specimen,  and  into  a  transmission  bar.  The  goal  of  this  test  is  to  subject  the  specimen  to  uniform 
uniaxial  stress  and  uniform  (biaxial)  strain  at  a  prescribed  axial  strain-rate.  These  uniform  conditions  can 
often  be  achieved  after  an  initial  “ring-up”  period  involving  multiple  wave  reflections  from  the  specimen- 
transmission  bar  and  specimen-incident  bar  interfaces.  Once  uniform  conditions  have  been  achieved  in 
the  specimen,  the  axial  stress,  axial  strain,  and  axial  strain-rate  histories  in  the  specimen  can  be  deduced 
from  strain  gage  measurements  on  the  elastic  bars:  the  axial  stress  is  proportional  to  the  strain  in  the 
transmission  bar,  and  the  axial  strain-rate  is  proportional  to  the  strain  in  the  incident  bar  generated  by  the 
wave  reflected  from  the  specimen  and  transmission  bar.  A  similar  approach  utilizes  a  torsional  wave  in 
the  Kolsky  bar  instead  of  a  compressive  wave  [Gilat  2000;  Hartley  et  al.  1985].  Typically,  a  number  of 
test  arc  performed  at  different  strain-rates,  and  the  data  are  presented  as  a  family  of  stress-strain  curves 
with  each  curve  corresponding  to  a  different  (approximately  constant)  strain-rate. 

The  standard  Kolsky  bar  techniques  work  well  on  metals,  composites,  and  stiff  polymers  and  have 
provided  useful  data  for  the  calibration  of  viscoplastic  constitutive  models  and  rate-dependent  failure 
models.  More  recently,  Kolsky  bar  tests  have  been  performed  on  gelatins  and  soft  biological  specimens 
with  the  goal  of  providing  data  for  the  calibration  of  viscoelastic  constitutive  models  for  these  mate¬ 
rials.  However,  considerable  difficulties  have  been  encountered  with  these  softer  specimens.  Radial 
accelerations  in  compression  tests  may  result  in  stress  states  that  are  not  uniaxial,  that  is,  radial  and 
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Figure  1.  Schematic  (not  to  scale)  of  the  experimental  setup  for  the  modified  torsional 
Kolsky  bar  described  in  [Nie  et  al.  2013].  Only  the  specimen  is  included  in  the  simulations. 

hoop  stresses  that  are  a  significant  fraction  of  the  axial  stress  but  cannot  be  measured  [Scheidler  and 
Kraft  2010;  Scheidler  et  al.  2011;  Warren  and  Forrestal  2010].  Studies  utilizing  compression  Kolsky  bar 
techniques  purportedly  removed  radial  inertia  effects  by  switching  from  disk-shaped  specimens  to  annular 
specimens  [Chen  and  Song  2011].  It  was  later  shown  that  annular  specimens  are  also  plagued  by  radial 
inertia  effects  and  that  the  observed  increases  in  axial  stress  at  high  strain-rates  were  falsely  attributed  to 
viscoelastic  properties  of  the  specimen  [Scheidler  and  Kraft  2010;  Sanborn  2011;  Sanborn  et  al.  2012]. 

Recently,  a  variation  of  the  torsional  Kolsky  bar  test  has  been  developed  for  studying  the  strain-rate 
sensitivity  of  very  soft  materials  and  applied  to  bovine  brain  tissue  [Nie  et  al.  2013;  Sanborn  et  al.  2012]. 
In  this  modified  torsion  test,  the  transmission  bar  is  replaced  with  a  torque  load  cell  in  order  to  generate  a 
stronger  signal  than  could  be  obtained  from  strain  gage  measurements  on  a  transmission  bar;  see  Figure  1 . 
The  faces  of  thin  annular  specimens  are  glued  to  the  incident  bar  and  the  load  cell.  Strain  gage  measure¬ 
ments  in  the  incident  bar  are  used  to  infer  the  rate  of  the  angle  of  twist  at  the  specimen-incident  bar  inter¬ 
face  and  consequently  a  volumetric  average  of  the  shear  strain-rate  in  the  specimen.  Integration  yields  the 
angle  of  twist  at  the  interface  and  consequently  a  volumetric  average  of  the  shear  strain  in  the  specimen. 
The  radially  averaged  shear  stress  at  the  specimen-load  cell  interface  is  deduced  from  the  torque  measured 
at  the  load  cell.  The  stress-strain  curves  in  [Nie  et  al.  2013]  are  obtained  by  plotting  this  radially  averaged 
shear  stress  at  the  load  cell  versus  the  volume  averaged  shear  strain  in  the  specimen.  Comparison  of 
stress-strain  curves  at  different  strain-rates  could  (in  principle)  yield  information  on  the  rate-sensitivity  of 
the  shear  stress.  Nie  et  al.  [2013]  performed  Kolsky  bar  tests  on  bovine  brain  tissue  at  an  average  shear 
strain-rate  of  700/s.  Analogous  quasistatic  tests  at  strain-rates  of  0.0 1/s  and  1/s  were  performed  for 
comparison,  and  the  results  appeared  to  indicate  a  substantial  increase  in  shear  stress  in  the  dynamic  test. 

The  effects  of  radial  inertia  in  this  torsional  Kolsky  bar  test  are  expected  to  be  much  lower  than  in 
a  compression  test  at  the  same  strain-rate.  However,  for  specimens  as  soft  as  brain  tissue,  the  effects  of 
axial  inertia  (i.e.,  torsional  wave  propagation)  may  be  significant.  The  shear  wave  speeds  in  the  softer 
specimens  are  so  low  that  even  for  thin  (i.e.,  1-2  mm  thick)  specimens  there  may  be  insufficient  time  for 
the  specimen  to  “ring-up”  to  an  axially  uniform  state  of  stress,  strain,  and  strain-rate.  In  particular,  if  the 
shear  strain  has  large  axial  variations,  then  the  shear  strain  at  the  specimen-load  cell  interface  may  not 
be  close  to  the  average  shear  strain  used  to  generate  the  stress-strain  curves  in  [Nie  et  al.  2013],  in  which 
case  their  “stress-strain  curves”  would  not  provide  valid  constitutive  data.  Likewise,  the  shear  strain-rate 
at  the  specimen-load  cell  interface  may  not  be  close  to  the  reported  average  shear  strain-rate. 

Torsion  and  compression  tests  on  soft  materials  differ  in  another  significant  way.  As  noted  above,  for 
compression  tests,  the  desired  uniaxial  stress  state  may  be  impossible  to  achieve  dynamically,  but  uniaxial 
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stress  can  always  be  achieved  quasistatically  if  frictional  effects  at  the  specimen  faces  are  negligible. 
However,  even  in  a  quasistatic  torsion  test,  the  desired  state  of  pure  shear  stress  is  impossible  to  achieve 
except  for  infinitesimal  deformations.  Significant  normal  stresses  can  develop  in  finite  deformation  qua¬ 
sistatic  torsion.  While  radial  stresses  can  be  reduced  by  the  use  of  annular  specimens,  axial  and  hoop 
stresses  may  not  be  negligible.  Furthermore,  the  glued  boundary  conditions  for  annular  specimens  can 
be  expected  to  result  in  axial  variations  in  the  stress  and  strain  state  even  in  a  quasistatic  test  since  radial 
deformations  arc  prohibited  at  the  ends  of  the  specimen  but  not  in  the  interior.  Thus,  unlike  quasistatic 
compression  tests,  the  stress  state  in  quasistatic  torsion  tests  on  soft  materials  is  already  quite  complicated, 
and  inertial  effects  in  the  corresponding  dynamic  test  will  further  complicate  these  stress  states.  Even 
if  the  full  stress  state  could  be  measured  in  a  dynamic  torsion  test,  it  would  be  difficult  to  sort  out  the 
role  of  inertial  effects  from  nonlinear  elastic  or  viscoelastic  effects  without  results  for  the  corresponding 
quasistatic  case  for  comparison.  This  problem  is  compounded  by  the  fact  that  in  the  torsional  Kolsky  bar 
test  only  the  torsional  (Oz)  stress  component  can  be  inferred  from  the  test  data,  and  this  is  known  only 
on  the  face  of  the  specimen  in  contact  with  the  load  cell. 

In  view  of  these  issues,  we  have  undertaken  numerical  simulations  of  both  dynamic  and  quasistatic 
torsion  tests  on  soft  and  nearly  incompressible  materials.  To  simplify  the  dynamic  simulations,  the 
incident  bar  and  the  load  cell  are  not  included;  they  are  replaced  with  appropriate  boundary  conditions 
at  the  specimen  ends.  Unlike  a  real  torsion  test,  the  simulations  provide  the  complete  stress  and  strain 
states  at  each  point  in  the  specimen  although  a  constitutive  model  must  necessarily  be  assumed.  A 
linear  elastic  constitutive  model  was  used  in  our  one-dimensional  simulations,  and  a  nonlinear  elastic 
constitutive  model  was  used  in  most  of  our  three-dimensional  simulations. 

Of  course,  soft  tissues  are  viscoelastic  (and  poroelastic  as  well),  and  the  goal  of  Kolsky  bar  tests 
on  these  inelastic  materials  is  to  quantify  their  strain-rate  sensitivity.  However,  performing  numerical 
simulations  with  elastic  constitutive  models  is  a  particularly  useful  way  of  revealing  inertial  effects  in 
dynamic  tests.  The  stress  in  an  elastic  material  is  rate-independent.  Thus,  any  differences  between 
quasistatic  and  dynamic  simulations  of  torsion  on  elastic  specimens  must  necessarily  be  due  to  inertial 
effects :  gradients  in  stress  and  strain  resulting  from  acceleration  of  the  specimen.  As  discussed  above, 
axially  nonuniform  shear  strains  in  the  specimen  would  invalidate  the  analysis  used  to  generate  the 
stress-strain  curves  from  the  experimental  data.  Use  of  viscoelastic  models  in  the  simulations  will  add 
strain-rate  effects  but  not  alter  the  presence  of  inertial  effects;  as  we  will  demonstrate,  the  main  issue 
here  is  the  slow  shear  wave  speed  in  the  specimen. 

Finally,  we  emphasize  that  the  purpose  of  this  paper  is  not  to  propose  or  examine  particular  constitutive 
models  for  brain  tissue  or  other  soft  materials  but  rather  to  investigate  the  validity  of  a  particular  dynamic 
test  that  was  developed  in  an  attempt  to  characterize  the  high  strain-rate  response  of  brain  tissue. 

The  paper  is  organized  as  follows.  First,  in  Section  2,  we  provide  an  overview  of  the  loading  pulses,  ma¬ 
terial  properties,  and  simulation  parameters.  Next,  in  Section  3,  we  discuss  our  one-dimensional  analysis 
and  simulations  of  torsional  waves  in  a  linear  elastic  specimen.  In  Section  4,  we  present  the  results  of  our 
quasistatic  analysis  and  simulations  of  a  full  three-dimensional  isotropic  nonlinear  elastic  specimen.  Here 
we  used  a  compressible  version  of  the  well-known  Mooney-Rivlin  model  for  incompressible  materials 
with  bulk  and  shear  moduli  chosen  to  yield  a  nearly  incompressible  material.  Section  5  contains  the 
results  of  our  three-dimensional  dynamic  simulations  of  the  torsion  test,  most  of  which  also  used  the 
Mooney-Rivlin  model.  However,  Section  5B  contains  a  brief  summary  of  three-dimensional  dynamic 
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simulations  that  used  a  nonlinear  viscoelastic  model  for  the  specimen.  Finally,  Section  6  contains  our 
major  conclusions  and  their  implications  on  the  experimental  data  collected  from  torsional  Kolsky  bar 
tests.  These  should  be  relevant  to  other  dynamic  test  techniques  that  attempt  to  generate  a  (locally) 
simple  shear  such  as  a  double-lap  shear  test  [Saraf  et  al.  2007;  Trexler  et  al.  2011]  or  a  circular  shearing 
test  [Nie  et  al.  2011].  In  Appendix  A,  we  provide  some  additional  details  of  the  numerical  methods  that 
were  used.  Appendix  B  contains  the  mathematical  details  of  the  loading  pulses  used  in  these  simulations. 


2.  Overview  of  specimen  properties,  loading  pulses,  and  simulations 

We  tailored  our  simulations  to  be  relevant  to  the  experimental  loading  conditions  and  specimen  geome¬ 
tries  considered  in  [Nie  et  al.  2013].  In  Sections  2A-2C,  we  discuss  how  we  arrived  at  our  material 
parameters,  boundary  conditions,  and  loading  pulse  shapes.  We  also  provide  an  overview  of  the  set  of 
simulations  in  Sections  2D  and  2E. 


2A.  Material  properties.  From  the  0.01  /s  and  1  /s  quasistatic  tests  on  bovine  brain  tissue  in  [Nie  et  al. 
2013],  we  estimated  small  strain  shear  moduli  of  /x  =  4.5  and  9kPa,  respectively.1  In  most  of  this  paper, 
we  use  a  shear  modulus  /x  =  8  kPa  as  representative  of  brain  tissue.2  This  is  near  the  upper  limit  of  the 
quasistatic  shear  moduli  reported  for  brain  tissue  in  the  literature,  the  lowest  values  being  around  2  kPa 
[Donnelly  and  Medige  1997;  Hrapko  et  al.  2006;  Miller  and  Chinzei  1997;  Nicolle  et  al.  2004;  Saraf  et  al. 
2007;  Shuck  and  Advani  1972],  In  view  of  the  high  water  content  of  brain  tissue,  we  assumed  that  the 
density  of  the  brain  tissue  specimens  was  that  of  water  (1  g/cm3).  This  specimen  density  was  used  in  all 
simulations  although  the  shear  modulus  was  varied  considerably  (as  described  in  Sections  2D  and  2E). 


2B.  Specimen  geometry  and  boundary  conditions.  The  annular  specimen  was  defined  as  a  hollow  cylin¬ 
der  of  outer  radius  R0  =  9.5  mm,  inner  radius  A',  =  7.35  mm,  and  relatively  short  length  L  —  1.7  mm, 
giving  a  small  length-to-diameter  ratio  of  0.09. 3  Let  R  denote  the  arithmetic  average  radius  and  A  the 
corresponding  difference  in  radii: 


-  Ra  +  Rj 

A  =  - 


A  = 


Ra  -  Ri 


—  Rn  —  R  —  R  —  Ri 


(2-1) 


so  that  Rj  =  R  —  A  and  R„  =  R  +  A.  Table  1  summarizes  the  geometrical  parameters  for  all  of  the 
simulations  performed. 

In  the  Kolsky  bar  tests,  the  faces  of  the  specimen  were  glued  to  the  incident  bar  and  load  cell,  and 
for  soft  specimens,  the  twist  at  the  specimen-load  cell  interface  should  be  negligible  compared  to  the 
twist  at  specimen-incident  bar  interface.  Thus,  for  all  simulations,  one  face  of  the  annular  specimen  was 
subjected  to  a  rigid  rotation  while  the  other  face  was  fixed  so  that  material  points  on  either  face  could 
not  undergo  axial  or  radial  displacements. 

Relative  to  a  cylindrical  coordinate  system,  let  (r,  9,  z)  and  (R,  0,  Z)  denote  the  coordinates  of  a 
material  point  in  the  deformed  configuration  and  undeformed  reference  configuration,  respectively.  The 

1  In  view  of  the  difficulties  in  extracting  an  initial  slope  from  stress-strain  data,  these  estimates  are  reasonably  close  to  their 
estimates  of  5.9  and  10.4kPa. 

2  An  instantaneous  elastic  shear  modulus  of  40kPa  was  used  in  the  viscoelastic  simulations  discussed  in  Section  5B.  The 
motivation  for  this  choice  is  discussed  there  and  in  Section  3D. 

3  No  information  is  provided  in  [Nie  et  al.  2013]  as  to  the  variability  in  specimen  dimensions  in  the  five  samples  they  tested. 
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R0  (mm)  Rj  (mm)  R  (mm)  A  (mm)  A/R  L  (mm)  L/2R0 
9.5  7.35  8.425  1.075  0.13  1.7  0.09 

Table  1.  Summary  of  geometrical  parameters  for  all  simulations  conducted. 


end  of  the  specimen  at  Z  =  L  is  taken  to  be  the  fixed  face  and  thus  corresponds  to  the  specimen-load 
cell  interface.  The  other  end  (Z  =  0)  corresponds  to  the  specimen-incident  bar  interface;  it  is  rotated  by 
an  angle  —  4i(f),  where  T  >  0  is  measured  in  radians  and  t  denotes  time.  Thus, 

0\z=0  (/)  =  0-^(0,  e\z=L(t)  =  ®.  (2-2) 


The  average  angle  of  twist  per  unit  length  is  given  by 


e\z=L(t)-o\z=o(t)  *(*)  n 

y(t)  = - 1 - =  — —  >  0. 


(2-3) 


L  L 

A  measure  of  the  average  shear  strain-rate  in  the  specimen  often  used  in  the  experimental  literature 
on  torsional  Kolsky  bar  tests  is 


■  0\z=L(t) -0\z=o(t)  B  V(t)  B  ■  - 

Ys(0  = - ; - R  =  ~^R  =  f(t)R, 


(2-4) 


where  a  superposed  dot  denotes  the  time  derivative  and  R  is  the  average  radius  introduced  above  [Gilat 
2000;  Nie  et  al.  2013].  This  average  shear  strain-rate  is  typically  inferred  from  strain  gage  measurements 
on  the  bar(s)  and  then  integrated  to  give  an  average  shear  strain  in  the  specimen: 


Ys(t)  = 


0\z=L(t )  ~  fl|z=o(0 

L 


R  - 


4AO 

L 


R  -  \j/{t)R  >  0. 


(2-5) 


While  other  measures  of  average  shear  strain  arc  often  used  in  the  theoretical  literature  (see  Section  4B), 
we  will  use  ys(t)  and  its  rate  in  order  to  simplify  comparison  with  the  experimental  data  in  [Nie  et  al. 
2013].  For  the  dynamic  simulations  considered  here,  we  applied  an  angular  displacement  at  Z  =  0  that 
delivered  the  desired  average  strain-rate  history  ys(t)  for  the  particular  simulation.  In  view  of  (2-2)  and 
(2-5),  this  displacement  history  is  given  by 


L 


0|z=o(O- ©  =  -*(*)  =  —=  MO- 

R 


(2-6) 


2C.  Loading  pulse  shapes.  We  use  the  term  loading  pulse  to  refer  to  the  average  str  ain-rate  history 
imposed  on  the  specimen.  The  key  features  of  a  typical  loading  pulse4  are  summarized  in  panel  a  of 
Figure  2.  These  loading  pulses  arc  intended  to  mimic  those  generated  in  Kolsky  bar  experiments.  The 
pulse  ramps  up  smoothly  over  the  time  interval  [0,  tr]  to  a  maximum  average  strain-rate  ymax,  which  is 
held  constant  for  the  plateau  duration  tp.  We  assume  a  symmetric  unloading  of  the  pulse;  that  is,  the 
unloading  time  is  equal  to  the  rise  time  tr  so  that  the  total  duration  td  of  the  loading  pulse  is  given  by 


td  —  2tr-\-tp.  (2-7) 

The  average  strain  accumulated  during  the  constant  strain-rate  portion  is  tpymax.  For  the  loading  pulse 
shapes  considered  here,  the  strains  accumulated  during  the  initial  rise  and  during  the  unloading  are 

4  A  complete  mathematical  description  of  this  idealized  loading  pulse  is  included  in  Appendix  B. 
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Figure  2.  Key  features  of  a  typical  average  strain-rate  history  for  the  dynamic  simula¬ 
tions  (panel  a).  Average  strain-rate  histories  (panel  b)  and  strain  histories  (panel  c)  for 
Simulations  I-IX. 


each  \trymax.  Thus,  the  total  average  strain  is 


Tmax  —  2(2  tr  Xmax)  T  A  Xmax  —  ilr  T  t  />)Y it 


(2-8) 


Nie  et  al.  [2013]  reported  stress-strain  curves  from  quasistatic  and  torsional  Kolsky  bar  tests  on  five 
samples  of  bovine  brain  tissue.  The  Kolsky  bar  tests  in  [Nie  et  al.  2013]  are  described  as  being  at  a 
strain-rate  of  700/s  with  a  pulse  duration  of  600  /x s.  The  final  average  shear  strain  for  their  tests  was 
about  0.3.  No  additional  details  of  the  loading  pulse  are  provided  with  the  following  exception.  The 
reader  is  referred  to  their  Figure  4  for  “a  typical  set  of  Kolsky  torsion  bar  result”.  This  figure  has  plots  of 
output  voltage  versus  time  for  the  two  gages.  The  only  information  that  we  could  infer  from  this  figure 
is  that  the  rise  time  of  the  loading  pulse  shown  was  about  250  /xs.5 

An  average  strain-rate  history  ys(0  consistent  with  the  key  features  described  above  was  obtained 
by  choosing  the  imposed  maximum  (average)  strain-rate  ymax  to  be  700/s,  the  rise  time  tr  to  this  700/s 
plateau  to  be  250 /xs,  and  the  maximum  (average)  strain  ymax  to  be  0.3.  These  conditions,  together  with 
(2-7)  and  (2-8),  yield  a  plateau  duration  of  tp  =  179  /xs  and  a  loading  pulse  duration  of  t,i  =  679  /xs.  The 
curves  for  the  one-dimensional  Simulations  /-///  in  Figure  2  show  the  average  strain-rate  history  ys(t) 
(panel  b )  and  the  corresponding  average  strain  history  ys(t)  (panel  c)  obtained  from  this  construction.  The 
other  two  curves  in  these  panels  are  the  average  strain  and  strain-rate  histories  for  the  one-dimensional 
Simulations  V,  VII,  and  IX  and  Simulations  IV,  VI,  and  VIII. 


5  From  the  discussion  in  [Nie  et  al.  2013],  it  is  not  clear  whether  the  data  in  this  figure  corresponds  to  any  of  the  tests  in  that 
paper  or  whether  the  data  is  taken  from  an  earlier  study  by  Nie  et  al.  [2011]  that  is  cited  frequently  in  the  2013  paper.  The  earlier 
paper  describes  a  circular  shearing  test,  but  it  is  referred  to  as  a  torsion  test.  Results  from  that  paper  are  cited  as  providing 
evidence  that  the  specimens  were  “undergoing  uniform  deformation”  for  the  700/s  tests  in  [Nie  et  al.  2013],  a  claim  that  is 
contradicted  by  the  results  presented  here. 
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Simulation 


parameter 

■  unit 

i 

II 

III 

IV 

V 

VI 

VII 

VIII 

IX 

Xmax 

0.3 

— 

— 

— 

— 

— 

— 

— 

— 

Xmax 

1/s 

700 

700 

700 

200 

500 

200 

500 

200 

500 

tr 

/XS 

250 

— 

— 

— 

— 

— 

— 

— 

— 

tp 

/xs 

179 

179 

179 

1250 

350 

1250 

350 

1250 

350 

td 

/xs 

679 

679 

679 

1750 

850 

1750 

850 

1750 

850 

P 

g/cm3 

1 

— 

— 

— 

— 

— 

— 

— 

— 

[l 

kPa 

800 

80 

8 

8 

8 

80 

80 

800 

800 

C 

mm//xs 

0.0283 

0.0089 

0.0028 

0.0028 

0.0028 

0.0089 

0.0089 

0.0283 

0.0283 

^travel 

/xs 

60.1 

190.1 

601 

601 

601 

190.1 

190.1 

60.1 

60.1 

Table  2.  Summary  of  material  and  loading  pulse  parameters  for  nine  one-dimensional 
simulations.  Long  dashes  ( — )  indicate  parameters  that  were  not  changed. 

2D.  One-dimensional  torsion  simulations.  The  one-dimensional  linear  elastic  simulations  described  in 
Sections  3C,  3D,  and  3F  were  conducted  with  the  primary  focus  on  the  effects  of  the  shear  modulus 
and  the  loading  pulse  on  the  quality  of  the  experiment.  Table  2  summarizes  the  parameters  specific  to 
the  nine  one-dimensional  simulations  in  Sections  3C  and  3D.  Columns  1  and  2  list  the  parameters  and 
respective  units,  and  the  remaining  columns  give  the  values  of  the  parameter  used.  The  shear  modulus 
was  varied  by  two  orders  of  magnitude:  800,  80,  and  8  kPa,  the  latter  being  on  the  order  of  brain  tissue,  as 
discussed  in  Section  2A.  These  shear  moduli  are  used  in  Simulations  I,  II,  and  III ,  respectively,  together 
with  the  700/s  loading  pulse  discussed  above.  In  Section  3C,  the  spatial  and  temporal  variation  of  the 
actual  shear  strain  and  shear  strain-rate  in  the  specimen  arc  compared  for  these  three  simulations. 

In  Section  3D,  we  compare  the  “stress-strain  curves”  for  the  same  three  shear  moduli  and  three  maxi¬ 
mum  (average)  strain-rates  (200,  500,  and  700/s)  for  a  total  of  nine  simulations  ( Simulations  I-IX).  The 
200/s  and  500/s  loading  pulses  had  the  same  rise  time  (250  /xs)  and  maximum  strain  (0.3)  as  the  700/s 
case  but  different  plateau  durations  and  total  pulse  durations  (see  Table  2  and  Figure  2).  The  “stress- 
strain  curves”  presented  here  arc  obtained  by  plotting  the  average  shear  stress  at  the  load  cell  versus  the 
average  shear  strain  in  the  specimen  as  in  the  experimental  results  reported  by  [Nie  et  al.  2013].  The  errors 
introduced  by  this  procedure  are  illustrated  by  comparison  with  the  actual  (linear)  stress-strain  curves  that 
would  have  been  obtained  if  the  specimen  were  undergoing  a  uniform  (i.e.,  pure  torsional)  deformation. 

In  Section  3F,  we  explore  the  plausibility  of  improving  the  quality  of  experimental  results  by  changing 
the  parameters  of  the  loading  pulse,  specifically  the  rise  time  tr  and  plateau  duration  tp.  In  these  studies, 
tr  took  on  the  thirteen  values  125,  250,  500,  750,  1000,  1250,  1500,  1750,  2000,  2250,  2500,  2750, 
and  3000  /xs  and  tp  took  on  the  eleven  values  1,  100,  200,  300,  400,  500,  600,  700,  800,  900,  and  1000  /xs 
for  a  total  of  143  additional  one-dimensional  simulations.  These  simulations  were  conducted  on  a  soft 
material  with  the  same  geometry  and  material  properties  as  in  Simulation  III. 

Additional  details  about  the  numerical  methods  for  the  one-dimensional  simulations  are  included  in 
Appendix  A. 

2E.  Three-dimensional  torsion  simulations.  In  all  of  our  three-dimensional  simulations,  stress-free 
boundary  conditions  were  used  on  the  outer  surface  R  =  R0  and  the  inner  surface  R  =  R,  of  the  annular 
specimen.  With  the  exception  of  the  viscoelastic  results  in  Section  5B,  the  constitutive  relation  for  the 
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specimen  was  a  compressible  version  of  the  nonlinear  elastic  Mooney-Rivlin  model  (see  Section  4).  We 
used  a  bulk  modulus  of  2.3  GPa  for  all  simulations;  this  is  the  bulk  modulus  of  water  and  approximates 
the  bulk  modulus  of  many  gelatins  and  soft  tissues  (such  as  brain  tissue). 

For  the  quasistatic  simulations  in  Section  4F,  the  shear  modulus  was  varied  by  almost  three  orders 
of  magnitude:  p  =  800,  80,  8,  and  2kPa.  The  corresponding  ratios  of  bulk  to  shear  modulus  varied 
from  2875  to  1.15  x  106  so  that  in  all  cases  the  material  was  nearly  incompressible.  For  each  value  of 
the  shear  modulus,  simulations  were  performed  for  three  values  of  the  nondimensional  parameter  co  in 
the  Mooney-Rivlin  model:  0.3,  0.6,  and  1  (neo-Hookean  model). 

The  three-dimensional  dynamic  simulations  ( Simulations  X-XVIII)  with  an  elastic  specimen  are  dis¬ 
cussed  in  Section  5A.  The  material  parameters  for  these  simulations  were  varied  in  an  analogous  way 
to  the  three-dimensional  quasistatic  simulations  except  that  the  p  =  2  kPa  case  was  not  considered.  The 
loading  parameters  were  the  same  as  in  the  one -dimensional  Simulations  /-///.  Section  5B  contains  a  brief 
discussion  of  some  three-dimensional  dynamic  simulations  that  used  a  nonlinear  viscoelastic  model  for 
the  specimen. 

Additional  details  about  the  numerical  methods  for  both  the  quasistatic  and  dynamic  three-dimensional 
simulations  arc  included  in  Appendix  A. 

3.  One-dimensional  treatment  of  dynamic  torsion 

Unlike  traditional  Kolsky  bar  tests  where  the  force  or  torque  is  measured  at  both  ends  of  the  specimen, 
in  the  modified  torsional  Kolsky  bar  test  of  [Nie  et  al.  2013],  the  torque  is  measured  by  a  load  cell  at  one 
end  of  the  specimen  only,  and  a  (radially  averaged)  shear  stress  at  that  end  of  the  specimen  is  extracted 
from  this  measurement.  The  test  provides  no  means  for  determining  whether  an  axially  uniform  state 
of  shear  stress  has  been  achieved  in  the  specimen.  This  would  not  be  an  issue  if  the  shear  strain  could 
also  be  measured  at  the  specimen-load  cell  interface,  in  which  case  a  valid  stress-strain  curve  could  be 
obtained.  However,  it  is  not  possible  to  measure  the  shear  strain  at  any  specific  location  in  the  specimen. 
Only  an  average  measure  of  shear  strain  in  the  specimen,  ys(t)  in  (2-5),  can  be  extracted  from  the  test 
data,  and  this  represents  a  volumetric  average.6 

As  pointed  out  in  the  introduction,  the  stress-strain  curves  reported  by  [Nie  et  al.  2013]  were  obtained 
by  plotting  the  (radially  averaged)  shear  stress  at  the  load  cell  versus  the  average  shear  strain  ys(t)  in  the 
specimen.  Clearly,  this  is  meaningful  only  if  ys(t)  is  close  to  the  actual  shear  strain  at  the  specimen-load 
cell  interface.  Such  a  condition  can  be  expected  to  hold  over  some  time  interval  only  if  the  shear  strain 
in  the  specimen  is  axially  uniform  on  that  time  interval.  In  an  effort  to  determine  when  a  state  of  uniform 
shear  strain  can  be  achieved,  we  have  performed  one-dimensional  simulations  of  dynamic  torsion  on  lin¬ 
ear  elastic  specimens.  Although  the  goal  of  the  torsional  Kolsky  bar  tests  is  to  measure  strain-rate  effects 
due  to  viscoelastic  properties  of  the  material,  a  study  of  one-dimensional  linear  elastic  wave  propagation 
in  the  torsion  specimens  reveals  the  difficulty  in  achieving  uniform  conditions  in  soft  specimens. 

In  Section  3A,  we  discuss  the  one-dimensional  wave  equation  that  is  solved  in  our  simulations.  Section 
3B  contains  a  brief  discussion  of  dynamic  equilibrium.  An  overview  of  the  simulations  in  Sections  3C, 
3D,  and  3F  was  given  in  Section  2D.  In  Section  3E,  we  consider  some  ways  of  improving  the  quality  of 
the  torsion  experiments. 

6  In  this  regard,  see  also  (3-11)  below. 
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3  A.  The  wave  equation  for  dynamic  torsion.  To  investigate  the  possibility  of  axially  nonuniform  shear 
strains  in  dynamic  torsion  tests,  we  consider  the  torsional  deformation 


so  that 


r  =  R,  6  —  6{Z,  t)  —  0  +  5(Z,  t),  z  =  Z 


8(z ,  t)  =  6{z ,  t)  -  0 


(3-1) 

(3-2) 


is  the  angular  displacement  at  the  axial  location  z  at  time  t.  Note  that  for  the  deformation  (3-1)  there  are 
no  axial  or  radial  displacements  in  the  specimen.  The  boundary  conditions  (2-2)  imply  that  the  angular 
displacements  at  the  ends  of  the  specimen  are 

8(L,t)  =  0,  5(0, 0  =  0(0,  =  =  (3-3) 


Recall  that  i Js(t)  is  the  average  angle  of  twist  per  unit  length  (see  (2-3)). 

We  assume  that  the  displacement  vector  u  and  its  spatial  gradient  V«  are  infinitesimal.  Let  e  denote 
the  infinitesimal  strain  tensor:  e  =  1(Vm  +  (V«)r).  Then  for  the  deformation  (3-1),  the  only  nonzero 
components7  of  u  and  e  relative  to  the  cylindrical  coordinate  system  are 


lie  =  r8(z ,  t ), 


=  e6>z  - 


1  dug 

2  ~di 


1  95 
—r  — . 

2  9  z 


We  also  assume  that  the  specimen  satisfies  the  isotropic  linear  elastic  constitutive  relation 


(3-4) 


T  =  A(tr  e)I  +  2  fie, 


(3-5) 


where  T  is  the  Cauchy  stress  tensor  and  the  constants  A.  and  fi  are  the  Lame  modulus  and  shear  modulus, 
respectively.  Then  for  the  deformation  (3-1),  the  only  nonzero  components  of  T  are 

dug  95 

Tze  =  Tgz  =  2fieez  =  fi—  =  fir  —  .  (3-6) 

dz  dz 


In  this  case,  the  equations  for  momentum  balance  in  cylindrical  coordinates  reduce  to 


dTgz  _  d2Ug 

dz  P  dt2 


(3-7) 


where  p  is  the  specimen  density.  Then  by  (3-6),  it  follows  that  ug  satisfies  the  one-dimensional  wave 
equation;  hence,  by  (3-4) i  and  (3-2),  5  and  6  satisfy  the  same  equation.8  For  example, 


925  _  2  925 
dt 2  dz2’ 


(3-8) 


That  is,  a  small  angular  displacement  will  propagate  axially  along  the  specimen  at  the  linear  elastic  shear 
wave  speed  c.9 


7  Here  and  elsewhere  in  the  paper,  we  use  the  physical  components  of  vectors  and  tensors,  that  is,  their  components  relative 
to  the  unit  basis  vectors  along  the  cylindrical  coordinate  directions.  See  [Graff  1975,  §A.9.1]  for  the  physical  components  of  e 
and  the  momentum  balance  equations  for  infinitesimal  strains. 

8  These  results  also  follow  front  strength-of-materials  approximations  [Graff  1975,  §2.6;  Wylie  1975,  §8.2]. 

9  A  shear  wave  propagating  into  an  initially  unstrained,  nonlinear  elastic  material  also  travels  at  the  linear  elastic  wave  speed 
[Truesdell  and  Noll  1965,  §73,  74,  77],  The  shear  wave  speed  may  change  after  the  wave  is  reflected  and  begins  to  travel  into 
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The  results  above  hold  for  both  disc-shaped  specimens  and  annular  specimens.10  In  either  case,  the 
outer  lateral  surface  r  =  R„,  and  for  annular  specimens  the  inner  lateral  surface  r  =  Rj,  is  stress-free 
since  Trr  =  TrH  =  Trz  =  0.  The  shear  strain  egz  and  the  shear  stress  Toz  vary  linearly  with  the  radius  as 
shown  above.  Restricting  attention  to  annular  specimens,  we  see  from  (3-6)  that  the  arithmetic  average 
of  the  shear  stress  at  the  outer  and  inner  lateral  surfaces  is 


where 


Taz\, — ^  +  TaAr-R.  —d8 

cr(z,  t)  =  9J’-R°  eJ,~R-  =  p R  =  py(z,  /), 

2  3  z 


-ds  _3  e 

y(z,  t )  =  R—  =  R  —  . 

az  az 


(3-9) 

(3-10) 


Note  that,  by  (3-4),  y(z,  t)  is  twice  the  arithmetic  average  of  f.uAz.  t)  at  the  outer  and  inner  lateral 
surfaces.11  In  interpreting  the  results  of  our  one-dimensional  simulations,  we  have  chosen  to  use  y(z,  t) 
as  a  measure  of  the  radially  averaged  shear  strain  at  the  axial  location  z  for  consistency  with  the  average 
shear  strain  ys(t)  used  in  the  experimental  literature.  Indeed,  by  (3-10),  (3-3),  and  (2-5),  ys(t)  is  the  axial 
average  value  of  y(z,t): 


1 

L 


y(z ,  t)dz  = 


-3.5 
R— dz 

3  z 


=  Ys(t). 


*[8(L,t)-8( 0,  t)]  =  Mt)R 


(3-11) 


Note  that  since  y(z,  t)  represents  a  radial  average  of  the  shear  strain  at  z,  ys(t)  represents  a  volumetric 
average  of  the  shear  strain  in  the  specimen. 

Recall  that,  in  the  dynamic  torsion  tests,  the  angle  of  twist  TO)  imparted  to  the  face  of  the  specimen 
(at  z  =  0)  is  typically  determined  in  such  a  way  that  a  desired  average  strain-rate  history  ys(t)  is  delivered. 
By  (3-2)  and  (2-6),  the  angular  displacement  history  necessary  to  deliver  that  average  strain-rate  history  is 

8(0,t)  =  -I=ys(t).  (3-12) 

R 

The  average  strain  histories  shown  in  Figure  2c,  when  substituted  into  (3-12),  yield  the  angular  displace¬ 
ment  boundary  conditions  that  are  actually  applied  at  the  loading  end  of  the  specimen  (z  =  0)  in  our 
simulations.  These  simulations  involve  the  numerical  solution12  of  the  one-dimensional  wave  equation 
(3-8)  for  8  with  zero  initial  conditions  and  boundary  conditions  (3-12)  and  (3-3) i  (i.e.,  8(L,  t)  =  0). 

A  pure  torsional  deformation  satisfying  the  boundary  conditions  (2-2)  is  given  by  (3-1)  with  angular 
displacement 

8(z,t)  =  f(t)(z  -  L),  f  >  0.  (3-13) 


a  strained  portion  of  the  material.  For  the  strains  considered  here,  nonlinearity  could  introduce  small  variations  in  wave  speeds 
that  would  be  minor  compared  with  other  effects  discussed  later. 

10  Because  of  the  small  length-to-diameter  ratio  of  the  specimens  considered  here,  namely  0.09,  the  more  common  terms 
“rod”  and  “tube”  do  not  seem  appropriate. 

1 1  For  the  annular  specimens  considered  here,  the  shear  strain  should  be  approximately  radially  uniform.  See  the  discussion 
of  pure  torsional  deformation  in  Section  4B. 

12  Additional  details  about  the  numerical  methods  for  the  one-dimensional  simulations  are  included  in  Appendix  A. 
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In  this  case,  by  (3-4),  the  only  nonzero  components  of  the  infinitesimal  strain  tensor  are 

cze  =  €ez  =  >  0.  (3-14) 

And  by  (3-10)  and  (2-5),  the  shear  strain  y  is  axially  uniform  and  equal  to  the  average  shear  strain  ys: 

y(z,  t)  —  \js(t)R  —  Ysijt)  for  pure  torsion.  (3-15) 

However,  a  state  of  pure  torsion  might  not  be  achieved  in  a  given  test  so  that  y(z,  t )  f  ys(t)  in  general. 
The  degree  to  which  y(z,t)  differs  from  its  measured  average  value  ys(t)  is  a  critical  point  of  this  paper 
and  is  discussed  in  the  following  sections. 


3B.  Dynamic  equilibrium.  In  the  experimental  literature,  the  terms  “force  equilibrium”,  “stress  equi¬ 
librium”,  “dynamic  stress  equilibrium”,  or  simply  “dynamic  equilibrium”  are  used  to  indicate  that  the 
axial  forces  (compression  tests)  or  torques  (torsion  tests)  on  the  two  ends  of  the  specimen  arc  balanced 
[Gray  2000;  Chen  and  Song  2011].  This  condition  is  clearly  necessary  (but  not  sufficient)  for  an  axially 
uniform  state  of  stress  and  strain  within  the  specimen,  and  in  the  experimental  literature,  these  uniform 
conditions  arc  often  (implicitly)  assumed  to  hold  on  any  time  interval  for  which  dynamic  equilibrium  is 
achieved.  We  will  use  the  term  “dynamic  equilibrium”  in  this  stronger  sense  of  axially  uniform  stress 
and  strain  states. 

Whether  or  not  dynamic  equilibrium  is  reached  depends  on  many  factors,  including  the  properties  of 
the  loading  pulse,  the  specimen  wave  speed,  and  the  specimen  length.  In  particular,  a  critical  factor  in 
achieving  dynamic  equilibrium  in  torsion  tests  is  the  travel  time  of  the  shear  wave: 


Gavel  —  —  L 

C 


(3-16) 


that  is,  the  time  it  takes  for  the  shear  wave  to  travel  from  one  end  of  the  specimen  to  the  other.  For  a 
uniform  strain  state  to  exist  in  the  material,  key  features  of  the  loading  pulse  such  as  the  rise  time  tr  and 
the  plateau  duration  tp  need  to  be  much  longer  than  the  travel  time  to  allow  the  specimen  to  “ring-up” 
to  dynamic  equilibrium  after  multiple  reflections. 

Table  2  lists  the  three  different  travel  times  for  Simulations  I— IX;  recall  that  we  assume  the  density  p 
is  that  of  water  and  the  specimen  length  L  is  1.7  mm.  For  the  values  of  the  shear  modulus  used  in  our 
dynamic  simulations  (8-800  kPa),  the  travel  time  ranges  from  601  /rs  to  60.1  /xs.  Since  the  travel  time 
is  the  delay  in  microseconds  between  the  initial  loading  of  the  specimen  at  z  =  0  and  the  arrival  of  the 
wave  at  the  load  cell  at  z  =  L,  the  way  this  travel  time  is  incorporated  into  the  analysis  of  experimental 
data  is  extremely  important  when  generating  stress-strain  curves  as  we  will  demonstrate  in  Section  3D. 
Assuming  that  the  shear  stress  o(L,  t)  at  the  load  cell  is  being  plotted  against  the  measured  averaged 
shear  strain  yft),  these  signals  would  have  to  be  shifted  in  time  to  synchronize  them.  However,  we  will 
show  that  this  procedure  can  introduce  considerable  errors. 


3C.  Failure  of  soft  specimens  to  reach  equilibrium.  In  this  section,  we  examine  the  spatial  and  temporal 
variation  of  the  strain  and  strain-rate  in  the  specimen  for  Simulations  I-III.  These  simulations  differ  only 
in  the  value  of  the  shear  modulus  p  used  for  the  specimen  (see  Table  2).  The  results  arc  summarized  in 
Figure  3.  False  color  images  of  the  strain-rate  (left  panels)  and  strain  (right  panels)  as  a  function  of  axial 
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Figure  3.  The  average  strain-rate  history  ys(t)  (panel  a)  and  corresponding  average 
strain  history  ys(t)  (panel  b)  imposed  on  the  specimen  for  Simulations  /-///.  The  first 
black  vertical  dashed  line  corresponds  to  the  rise  time  tr  to  constant  strain-rate.  The 
second  black  vertical  dashed  line  is  the  time  (tr  +  tp)  at  which  the  strain-rate  starts  to 
ramp  down  after  a  constant  strain-rate  plateau  of  duration  tp.  The  remaining  panel  pairs 
(c-d,  e-f,  and  g-h)  are  false  color  images  of  the  corresponding  actual  strain-rates  y(z,  t) 
and  strains  y{z,  t)  in  specimens  with  shear  moduli  /x  =  800,  80,  and  8kPa,  respectively. 
Note  that  the  color  scale  varies  for  the  three  groups  of  plots.  The  vertical  magenta  lines 
correspond  to  ftraVeb  which  depends  on  the  shear  modulus. 


position  (vertical  axis)  and  time  are  presented  in  panels  c-h.  Panels  a  and  b  show  the  average  strain-rate 
and  strain  histories  that  are  imposed  on  the  specimen. 

Figure  3  illustrates  the  degree  to  which  dynamic  equilibrium  has  (or  has  not)  been  achieved.  Shear 
waves  are  initiated  at  the  specimen-bar  interface  (z  =  0)  and  propagate  to  the  Load  Cell  (z  =  L).  In  each 
panel,  the  first  vertical  black  dashed  line  indicates  the  rise  time  tr.  The  two  vertical  black  dashed  lines 
taken  together  represent  the  window  of  time  during  which  the  imposed  average  strain-rate  is  constant. 
The  vertical  magenta  dashed  line  indicates  the  travel  time  and  thus  the  arrival  of  the  shear  wave  at  the 
Load  Cell,  which  increases  with  decreasing  shear  modulus  (see  (3-16)).  Thus,  the  black  dashed  lines 
indicate  the  key  time  scales  of  the  loading  pulse  and  the  magenta  line  represents  the  critical  travel  time 
of  the  specimen.  Ideally,  the  magenta  line  should  lie  well  before  both  black  lines. 

The  results  for  /x  =  800  kPa  ( Simulation  7)  are  shown  in  panels  c  and  d  of  Figure  3.  These  results 
are  close  to  what  one  might  expect  for  a  valid  torsional  Kolsky  bar  experiment.  The  strain-rate  (panel  c) 


NONUNIFORM  SHEAR  STRAINS  IN  TORSIONAL  KOLSKY  BAR  TESTS  ON  SOFT  SPECIMENS 


527 


increases  to  an  approximately  axially  uniform  (intended)  plateau  value  of  700/s  before  dropping  back 
down  to  a  zero  strain-rate  but  nonzero  final  strain  (panel  d).  From  the  false-color  representation,  one  can 
conclude  that  both  the  strain  and  strain-rate  are  fairly  uniform  along  the  length  of  the  specimen  at  any 
instant.  This  is  indicated  by  a  nearly  uniform  color  in  the  vertical  direction  at  any  instant.  The  ringing 
that  can  be  seen,  especially  in  the  strain-rate,  is  one  of  the  fundamental  modes  of  the  specimen.  There  is 
moderate  variability  in  the  strain-rate  due  to  the  ringing,  but  the  intended  strain-rate  and  maximum  strain 
nearly  match  their  imposed  average  values  plotted  in  panels  a  and  b.  Thus,  panels  c  and  d  in  Figure  3  are  a 
representation  of  a  strain  state  where  dynamic  equilibrium  is  nearly  achieved.  If  dynamic  equilibrium  was 
actually  reached,  the  color  representation  of  y(z,  t)  and  y(z,t)  would  not  change  in  the  vertical  direction. 

The  results  for  the  smaller  shear  modulus,  /x  =  8kPa  ( Simulation  III),  are  shown  in  panels  g  and  h 
of  Figure  3.  Clearly,  the  strain  and  strain-rate  in  Simulations  III  and  I  arc  considerably  different  both  in 
magnitude  and  uniformity.  Recall  that  a  shear  modulus  of  8  kPa  is  on  the  order  of  brain  tissue,  so  the 
specimen  stiffness  in  this  case  most  closely  resembles  the  stiffness  of  the  brain  tissue  specimens  in  the 
tests  of  [Nie  et  al.  2013].  Since  the  linear  elastic  shear  wave  speed  is  proportional  to  the  square  root  of  the 
shear  modulus  and  the  shear  modulus  in  Simulation  III  is  one  hundred  times  smaller  than  in  Simulation  I, 
the  wave  speed  is  ten  times  slower  than  in  Simulation  I.  The  travel  time  is  therefore  ten  times  longer, 
i.e.,  601  /xs  instead  of  60  /xs  (compare  the  magenta  lines  for  panels  c  and  d  against  panels  g  and  h  in 
Figure  3).  For  the  /x  =  8  kPa  case,  the  travel  time  is  so  close  to  the  pulse  duration  ( tj  =  679  /xs)  that  the 
shear  strain  and  shear  strain-rate  within  the  specimen  are  nonuniform  traveling  pulses.  Thus,  the  same 
loading  conditions  that  produced  axially  uniform  shear  strain  and  shear  strain-rate  in  the  /x  =  800  kPa 
case  resulted  in  highly  nonuniform  shear  strain  and  shear  strain-rate  within  the  8  kPa  specimen. 

For  the  /x  =  8  kPa  case,  it  is  grossly  inaccurate  to  describe  the  strain  and  strain-rate  at  any  location  in 
the  specimen  by  their  average  values  ys  and  ys  although  these  averages  arc  all  that  can  be  inferred  from 
test  measurements.  The  following  features  illustrate  the  degree  of  nonuniformity  in  the  specimen: 

(1)  During  the  first  600  /xs  of  the  simulation,  that  is,  prior  to  reflection  of  the  shear  wave  from  the  Load 
Cell,  the  peak  strain-rate  and  peak  strain  in  the  specimen  arc  about  3200/s  and  0.42,  respectively. 
These  peak  values  occur  at  different  axial  locations  at  different  times  (see  panels  g  and  h)  and  well 
exceed  the  maximum  average  strain-rate  (ymax  =  700/s)  and  maximum  average  strain  (ymax  =  0.3). 

(2)  The  largest  strain  achieved  in  the  specimen  during  the  simulation  is  about  0.84.  This  occurs  near 
the  end  Z  =  1 .7  mm  (i.e.,  adjacent  to  the  Load  Cell)  at  about  900  /xs  (panel  h).  This  peak  strain  is 
nearly  three  times  larger  than  the  average  strain  ys  at  900  /xs,  which  has  already  reached  its  final 
value  of  ymax  =  0.3. 

(3)  The  largest  strain-rate  achieved  in  the  specimen  during  the  simulation  is  about  5800/s.  This  also 
occurs  at  the  Load  Cell  end  of  the  specimen  (Z  =  1.7  mm)  at  about  730  /xs,  by  which  time  the 
average  strain-rate  ys  has  unloaded  to  zero.  This  peak  strain-rate  is  more  than  eight  times  higher 
than  the  maximum  average  strain-rate  (700/s)  imposed  on  the  specimen. 

(4)  Large  positive  strain-rates  and  a  zero  average  strain-rate  at  the  same  instant  imply  negative  strain- 
rates  of  large  absolute  value  must  also  be  experienced  in  some  portions  of  the  specimen.  This  is 
evidenced  by  the  blue  bands  in  panel  g. 

(5)  The  doubling  of  the  peak  strains  (from  0.42  to  0.84)  and  the  near  doubling  of  the  peak  strain-rates 
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(from  3200/s  to  5800/s)  is  due  to  the  reflection  of  the  shear  wave  from  the  fixed  end  (Z  =  1.7  mm) 
beginning  at  601  /xs.13 

(6)  As  an  example  of  the  large  axial  gradients,  the  strain  at  300 /xs  decreases  from  0.275  to  0  over  a 
distance  of  0.3  mm.  However,  the  change  in  the  angle  of  twist  over  this  0.3  mm  length  is  quite 
small:  roughly  0.005  radians  or  about  0.3°.  This  is  well  below  the  resolution  of  the  high  speed 
camera  images  of  the  specimen  in  [Nie  et  al.  2011]  that  were  often  cited  in  [Nie  et  al.  2013]  as 
evidence  that  the  specimens  undergo  uniform  deformation. 

The  results  for  /x  =  80kPa  ( Simulation  /)  are  shown  in  panels  e  and  /  of  Figure  3.  Even  for  this 
intermediate  case,  dynamic  equilibrium  is  not  reached.  Further  evidence  of  this  fact  is  given  by  the  stress- 
strain  curves  for  /x  =  80  kPa  in  the  next  section.  Clearly,  achieving  dynamic  equilibrium  in  specimens 
with  shear  moduli  substantially  less  than  800  kPa  can  be  quite  problematic. 

3D,  “Stress-strain  curves”  for  linear  elastic  specimens.  Up  to  this  point,  we  have  used  the  linear  elastic 
simulations  only  for  predictions  of  the  strain  and  strain-rate  within  the  specimen.  In  this  section,  we 
use  the  results  of  those  simulations  to  construct  “stress-strain  curves”  at  several  “constant  strain-rates”, 
analogous  to  the  way  that  test  data  is  typically  analyzed  and  presented  in  the  experimental  literature. 

Let  OLoadCeiiCO  a°d  /Load CeiiO  denote  the  (radially  averaged)  shear  stress  and  shear  strain  at  the 
Load  Cell : 

OLoadCell  (t)  =  cr(L,t),  Hoad  CellO)  =  Y  (L  -  0 .  (3-17) 

In  the  torsional  Kolsky  bar  test  of  [Nie  et  al.  2013],  the  torque  is  measured  at  the  load  cell  and  crLoadCeii  (7 ) 
can  be  extracted  from  this  measurement.  However,  HoadCeiiU)  cannot  be  measured  in  this  or  any  other 
torsional  Kolsky  bar  test  —  only  the  average  shear  strain  ys(t)  can  be  inferred  as  discussed  in  Section  2B. 
Consequently,  the  stress-strain  curves  in  [Nie  et  al.  2013]  are  plots  of  OLoadCell  versus  ys.  But  while 
ys(t)  starts  increasing  at  t  —  0,  that  is,  as  soon  as  the  twist  is  applied  to  the  end  z  =  0,  there  is  a  delay 
(equal  to  the  travel  time)  before  a  nonzero  strain  is  experienced  by  the  specimen  at  the  load  cell  (see 
Figure  3).  Since  this  travel  time  is  a  significant  portion  of  the  test  duration  for  the  soft  specimens  of 
interest,  it  might  seem  reasonable  to  time-shift  the  <T|  oad  CeiiO)  history  in  order  to  synchronize  it  with 
the  ys(t)  history  (or  vice  versa)  prior  to  plotting  one  against  the  other.  On  the  other  hand,  a  time  shift 
is  unnecessary  if  the  specimen  is  in  dynamic  equilibrium ,  a  condition  that  Nie  et  al.  [2013]  assumed  to 
hold  in  their  tests.  Since  the  issue  of  travel  time  in  the  specimen  was  not  discussed  in  their  paper,  it  was 
not  deal-  to  us  whether  a  time  shift  had  been  performed  in  generating  their  stress-strain  curves. 

In  this  section,  we  utilize  the  results  from  Simulations  I-IX  (see  Table  2)  to  generate  “stress-strain 
curves”  (i.e.,  OLoadCell  versus  ys)  with  and  without  time  shifts.  We  compare  the  results  of  both  approaches 
in  order  to  illustrate  the  different  apparent  “strain-rate  effects”  that  each  produces.  Recall  that  the  simu¬ 
lations  solve  the  one-dimensional  wave  equation  (3-8)  for  the  angular  displacement  8  with  the  boundary 
condition  at  z  =  0  determined  by  the  appropriate  average  strain  history  (see  (3-12)  and  Figure  2).  Six 
sets  of  “stress-strain  curves”  are  plotted  in  Figure  4.  Solid  colored  curves  arc  shear  stress  measured  at 
the  Load  Cell  plotted  against  the  average  shear  strain  ys.  The  colored  curves  in  the  left  column  do  not 
involve  any  time  shifts  of  the  data;  i.e.,  they  are  parametric  plots  of  OLoadCell (0  versus  Ks(0-  The  colored 

13  Another  doubling  would  be  expected  to  occur  at  the  specimen-bar  interface  (Z  =  0)  if  the  simulation  had  been  continued 
to  longer  times. 
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curves  in  the  right  column  are  parametric  plots  of  OLoad  Ceii(t  —  travel)  versus  ys(i);  i.e.,  the  Load  Cell  data 
has  been  time-shifted  to  synchronize  it  with  the  average  strain  data.  Each  row  of  plots  is  for  a  unique 
shear  modulus  (800,  80,  and  8kPa),  and  each  color  within  a  panel  represents  a  simulation  taken  at  a 
different  intended  maximum  average  strain-rate:  ymax  =  700/s  (blue),  500/s  (green),  and  200/s  (red). 
The  average  strain-rate  histories  are  shown  in  panel  b  of  Figure  2.  The  maximum  strain  ymax  and  rise 
time  were  held  constant  at  0.3  and  250  /xs,  respectively,  and  only  the  plateau  time  and  the  total  duration 
were  changed  to  accommodate  the  various  strain-rates  (see  Table  2). 

In  the  simulations  the  stress  and  strain  at  the  Load  Cell  are  linearly  related:  setting  z  =  L  in  (3-9)  gives 

OLoad  Cell  (0  =  PkLoad  Cell(0-  (3-18) 

Thus,  regardless  of  the  strain  and  strain-rate  histories  experienced  by  the  specimen,  a  parametric  plot 
of  OLoad  Cell  it )  versus  yi.oad  Ceil  (0  should  be  linear  with  slope  /x.  We  cannot  actually  plot  the  linear  relation 
in  (3-18)  in  Figure  4  since  the  variable  on  the  horizontal  axis  is  ys(t),  which  (as  shown  in  Section  3C) 
is  generally  not  equal  to  /Load  Ceii(0-  However,  if  the  shear  strain  were  uniform,  then  ys  would  coincide 
with  /Load  Cell ,  in  which  case  the  shear  stress  at  the  Load  Cell  would  vary  linearly  with  ys:  OLoad  Cell  =  /xys- 
This  lineal-  curve  is  shown  by  the  black  dashed  lines  in  each  panel  in  Figure  4.  It  represents  the  correct 
stress-strain  curve  against  which  the  colored  “rate-dependent”  curves  should  be  compared.  That  is,  while 
each  of  the  colored  curves  are  valid  plots  of  OLoad  Cell  versus  ys  (with  or  without  time  shifts),  they  represent 
the  crLoad Ceil  versus  yn>ad  Cell  curve  only  to  the  extent  to  which  they  agree  with  the  black  dashed  lines 
in  Figure  4.  Recall  that  the  shear  strain  y(z,  t)  is  uniform  for  a  specimen  undergoing  a  pure  torsional 
deformation  (see  the  discussion  at  the  end  of  Section  3A).  Hence,  the  black  dashed  lines  represent  the 
stress-strain  curves  that  would  be  obtained  for  pure  torsion. 

Panel  a  in  Figure  4  corresponds  to  a  shear  modulus  of  800  kPa  without  a  time  shift.  Note  that  all 
four  curves  lie  nearly  on  top  of  each  other.  Here  the  blue  curve  represents  the  results  from  Simulation  I. 
As  discussed  earlier,  this  was  considered  a  “valid”  Kolsky  bar  test,  and  that  conclusion  is  confirmed 
here.  Panel  h  is  the  same  set  of  stresses  shifted  by  the  travel  time  (60  /xs).  This  has  the  appearance  of 
valid  Kolsky  bar  data.  Note,  however,  that  all  these  curves  are  vertically  offset  from  the  true  (dashed) 
response,  implying  an  apparent  rate-dependence.  Despite  the  vertical  offset,  their  slopes  capture  the  true 
lineal'  behavior  of  the  material. 

We  emphasize  that  if  a  linear  elastic  specimen  is  in  dynamic  equilibrium,  there  should  be  absolutely 
no  strain-rate  effects.  However,  Figure  4  clearly  shows  that,  for  the  smaller  shear  moduli  of  80  and 
8  kPa,  the  apparent  “rate  effects”  are  substantial  regardless  of  whether  the  Load  Cell  data  is  time-shifted. 
It  is  unclear  what  information  could  have  been  inferred  from  these  plots  if  the  true  linear  stress-strain 
relation  were  unknown  with  the  exception  that,  for  /x  =  80  kPa,  a  linear  fit  to  the  unshifted  curves  or  to 
the  time-shifted  200/s  curve  might  return  a  slope  close  to  /x. 

When  the  shear  stress  is  plotted  without  a  time  shift,  there  is  an  initial  strain  interval  for  which  the 
stress  is  zero.  The  strain  duration  of  this  initially  flat  portion  of  the  stress-strain  curve  increases  with 
the  strain-rate  and  decreases  with  increasing  shear  modulus.  These  initially  flat  regions  are  present  even 
for  the  /x  =  800  kPa  case  (panel  a)  although  they  cannot  be  observed  on  the  scale  of  that  figure.  This 
nonphysical  behavior  is  most  profound  for  the  softest  specimen  (/x  =  8  kPa,  panel  e).  In  particular,  for 
the  700/s  case  in  that  panel,  there  is  no  measured  stress  until  ys(t)  is  close  to  0.3.  Recall  that  it  takes 
601  /xs  for  the  shear  wave  to  arrive  at  Load  Cell  whereas  the  average  strain  has  reached  its  final  value 
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Without  Time  Shift  With  Time  Shift 


Figure  4.  Panel  pairs  a-b,  c-d,  and  e-f  show  apparent  “strain-rate  effects”  in  three 
linear  elastic  materials  with  shear  moduli  /z  =  800,  80,  and  8  kPa,  respectively.  Note  the 
differences  in  the  vertical  scales  for  the  three  groups  of  panels.  Solid  colored  curves  are 
shear  stress  measured  at  the  Load  Cell  plotted  against  the  average  shear  strain  ys.  The 
black  dashed  lines  correspond  to  the  stress-strain  curve  that  results  for  a  pure  torsional 
deformation  in  dynamic  equilibrium. 

of  0.3  by  679  /zs  (the  pulse  duration).  For  t  >  ftrave]  =  601  /zs,  the  strain  at  the  Load  Cell  increases  to 
roughly  0.8  (see  panel  h  in  Figure  3)  and  subsequently  falls  to  zero  strain,  and  of  course,  the  stress  at 
the  Load  Cell  is  proportional  to  this  strain.  Thus,  a  parametric  plot  of  OLoadCeitCO  versus  ys(t)  becomes 
essentially  a  vertical  line  at  ys  =  0.3,  where  a  single  average  shear  strain  maps  to  infinitely  many  stresses. 

When  the  /z  =  80  and  8  kPa  data  is  plotted  with  a  time  shift  (panels  d  and /  in  Figure  4),  the  data  appears 
more  reasonable  but  is  still  completely  erroneous.  In  particular,  the  “initial”  shear  moduli  (that  is,  the 
slopes  of  the  stress-strain  curves  at  the  origin)  are  substantially  higher  than  the  actual  shear  modulus  of  the 
linear  elastic  specimen,  which  is  the  slope  of  the  dashed  curve.  For  strains  up  to  0.05,  our  700/s  curves 
(for  /z  =  8  and  80kPa)  agree  qualitatively,  though  not  quantitatively,14  with  the  700/s  stress-strain  curves 
for  bovine  brain  tissue  in  Figures  5  and  6  of  [Nie  et  al.  2013].  In  particular,  their  stress-strain  curves  have 
no  initially  flat  region  of  the  type  in  panels  c  and  e  in  our  Figure  4.  Based  on  these  observations,  it  would 
appear  that  Nie  et  al.  [2013]  time-shifted  the  Load  Cell  data  to  synchronize  it  with  the  average  strain 
data  in  order  to  generate  the  stress-strain  curves  for  the  “700/s”  tests  that  they  reported.  But  regardless 

14  Our  linear  elastic  simulations  for  /max  =  700/s  were  not  intended  or  expected  to  fit  the  700/s  stress-strain  curves  in  [Nie 
et  al.  2013],  Their  quasistatic  curves  show  large  nonlinear  elastic  effects  that  would  be  expected  to  persist  in  dynamic  tests. 
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of  whether  they  time-shifted  their  Load  Cell  data,  the  results  in  this  section  indicate  that  the  stress-strain 
curves  reported  for  bovine  brain  tissue  in  [Nie  et  al.  2013]  most  likely  have  large  systematic  errors  in  the 
strains.  The  results  in  the  next  paragraph  lend  further  support  for  this  conclusion. 

The  700/s  curve  in  Figure  6  of  [Nie  et  al.  2013]  is  an  average  of  the  700/s  stress-strain  curves  for 
their  five  samples.  From  this  curve,  they  estimated  an  (average)  initial  shear  modulus  of  53.5  kPa.  If  this 
were  an  accurate  estimate  of  the  instantaneous  elastic  shear  modulus  of  a  viscoelastic  specimen,  then 
the  shear  wave  speed  in  the  specimen  would  be  controlled  by  this  modulus,  giving  a  travel  time  in  the 
specimen  of  232  /zs.  Our  /z  =  80  kPa  elastic  simulations  are  closest  to  this  case  with  an  even  shorter 
travel  time  of  190  /zs.  However,  as  panels  c  and  d  in  Figure  4  and  panels  e  and  /  in  Figure  3  show, 
dynamic  equilibrium  is  not  reached  even  for  this  stiffer  case. 

It  is  clear  that  the  apparent  “rate  effects”  discussed  here  are  not  due  to  viscoelastic  or  viscoplastic 
behavior  in  the  specimen.  Rather,  they  are  due  to  the  fact  that  the  average  shear  strain,  which  is  used 
to  construct  the  stress-strain  curves,  is  significantly  different  from  the  strain  at  the  location  where  the 
stress  is  measured.  However,  in  a  Kolsky  bar  test  on  specimens  that  exhibit  rate-dependent  behavior, 
the  distinction  between  the  actual  material  response  and  apparent  “rate  effects”  would  be  lost.  Thus, 
one  might  falsely  attribute  the  apparent  rate  effects  discussed  above  to  either  viscoelastic  or  viscoplastic 
properties  of  the  specimen. 

After  this  portion  of  the  paper  had  been  completed,  private  communications  with  Xu  Nie  (the  first 
author  of  [Nie  et  al.  2013])  suggested  some  possibilities  that  were  not  considered  above.  Nie  stated 
that  in  generating  their  stress-strain  curves,  they  did  not  time-shift  the  Load  Cell  data  to  synchronize 
it  with  the  average  strain  data.15  Furthermore,  he  estimated  that  the  travel  time  in  the  specimen  was 
approximately  30  /zs.  This  estimate  was  obtained  from  the  difference  between  the  arrival  time  of  the 
loading  wave  at  the  specimen  and  the  onset  of  a  signal  at  the  Load  Cell.  These  times  were  estimated 
from  oscilloscope  traces  from  one  test.16  Clearly,  a  travel  time  of  30 /zs  is  inconsistent  with  the  results 
presented  here  for  a  specimen  as  soft  as  brain  tissue.  Indeed,  it  would  imply  an  instantaneous  elastic 
shear  modulus  of  3.2  MPa,  which  is  totally  unreasonable  for  brain  tissue. 

On  examining  a  snapshot  of  the  oscilloscope  traces  provided  to  us  by  Xu  Nie,  we  noticed  that  there 
was  a  discontinuity  in  the  slope  of  the  Load  Cell  signal  at  about  260/zs.17  Subsequently,  the  Load  Cell 
signal  exhibits  a  substantial  increase  in  noise.  We  suspect  that  the  travel  time  of  the  torsional  shear  wave 
in  the  specimen  was  (approximately)  260  /zs  and  that  the  onset  of  the  Load  Cell  signal  at  30  /zs  is  due  to 
some  other  source.  Note  that  a  travel  time  of  260  /zs  would  imply  an  instantaneous  elastic  shear  modulus 
of  42kPa,  which  is  not  unreasonable.  The  early  Load  Cell  signal  would  also  explain  why  the  stress-strain 
curves  in  [Nie  et  al.  2013]  do  not  exhibit  an  initially  flat  region  even  though  they  claim  that  the  data  was 
plotted  without  a  time  shift.  Indeed,  no  such  region  is  observable  in  panel  a  of  Figure  4  even  though  the 
travel  time  in  that  case  was  60  /zs.  If  our  conjecture  is  correct,  then  the  shear  stress  inferred  from  the 

15  They  did,  of  course,  perform  the  necessary  time  shift  of  the  strain  gage  data  to  account  for  the  time  it  takes  the  shear  wave 
to  propagate  from  the  strain  gage  to  the  specimen  end  of  the  incident  bar. 

16  Nie  points  out  (and  we  agree)  that  such  estimates  for  arrival  and  onset  times,  and  particularly  for  their  differences,  are 
subject  to  errors  that  are  difficult  to  bound.  Our  estimate  of  the  actual  travel  time  of  the  shear  wave,  and  consequently  our 
estimate  for  the  instantaneous  elastic  shear  modulus  for  the  specimen,  are  subject  to  even  more  uncertainty.  However,  it  is  only 
the  order  of  magnitude  of  these  estimates  that  are  essential  for  the  arguments  given  here. 

1 7  To  simplify  the  discussion,  we  take  the  arrival  time  of  the  loading  wave  at  the  specimen  to  be  r  =  0  as  has  been  assumed 
previously  this  paper. 
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Load  Cell  signal  would  be  erroneous  from  30-260  /zs  (and  hence  for  the  smaller  strains),  and  possibly 
for  later  times  as  well,  depending  on  how  long  this  spurious  signal  persisted. 

A  possible  explanation  for  the  early  onset  of  the  Load  Cell  signal  is  the  presence  of  bending  or 
longitudinal  waves  in  the  incident  bar,  which  would  travel  faster  than  the  torsional  shear  wave.  Another 
possible  explanation  is  the  longitudinal  precursor  wave  that  is  generated  in  the  specimen  by  the  shearing 
deformation  from  the  torsional  wave.  This  precursor  wave  is  a  real  effect  as  discussed  in  the  analysis  of 
our  three-dimensional  dynamic  simulations  in  Section  5.  But  it  is  not  clear  to  us  whether  any  of  these 
waves  could  even  generate  a  measurable  signal  in  the  torque  sensor,  so  at  the  moment,  the  early  onset  of 
this  signal  remains  unresolved. 

3E,  Increasing  the  quality  of  the  test.  In  the  previous  sections,  we  examined  the  axially  non  uniform 
strain  and  strain-rate  states  that  may  exist  in  the  specimen  in  torsional  Kolsky  bar  tests.  For  sufficiently 
soft  specimens,  large  systematic  error  in  the  stress-strain  curves  is  associated  with  the  loading  conditions 
and  the  specimen  lengths  considered.  In  this  section,  we  discuss  some  alterations  to  the  experiment  that 
might  alleviate  this  error.  Some  dimensional  analysis  is  particularly  useful  here  and  can  indicate  how  to 
improve  the  experimental  setup. 

Using  parameters  specific  to  our  loading  pulse,  by  (2-7)  and  (2-8),  we  see  that  the  total  average  strain 
and  the  pulse  duration  can  be  expressed  as 

Tmax  =  (tr  +  tp)ymSiX,  td  —  2tr  +  tp  —  tr  +  — - =  2- - tp.  (3-19) 

Tmax  Kmax 

If  we  take  ymax  as  fixed,  we  need  only  specify  two  of  the  three  parameters  tr,  tp,  and  ymax  in  order  to 
determine  the  four  parameters  tr,  tp,  ymax,  and  trj.  We  choose  to  track  changes  by  altering  the  rise  time  tr 
and  the  plateau  duration  tp. 

Two  quality  measures  of  the  test  conditions  arc  the  ratio  Qr  of  the  rise  time  to  the  travel  time  and  the 
ratio  Qp  of  the  plateau  time  to  the  travel  time: 

tr  Cty  tr  I 

flravel  L  Ly  p  * p 

Qr  can  be  thought  of  as  the  number  of  wave  reflections  that  occur  during  the  initial  rise  time,  and  Qp  can 
be  thought  of  as  the  number  of  reflections  that  occur  during  the  plateau  period.  A  high  quality  test  would 
be  one  for  which  both  Qr  dp  1  and  Qp  3>  1 .  We  emphasize  that  these  measures  are  dimensionless  and 
depend  on  the  details  of  the  loading  pulse,  the  specimen  length,  and  the  shear  wave  speed  of  the  specimen. 

For  a  given  range  of  values  of  the  shear  wave  speed  c,  the  parameters  L,  tj,  tr,  tp,  ymax,  or  ymax 
can  be  chosen  to  ensure  a  large  quality  factor.  The  //  =  800  kPa  case  in  Figure  4  (panel  a)  gave  the 
most  accurate  results:  here  Qr  =  4.2  for  all  three  cases,  and  Qp  =  2.98,  5.82,  and  20.8  for  ymax  =  700/s, 
500/s,  and  200/s,  respectively.  For  our  particular  examples,  transitioning  from  //  =  800  kPa  to  p  =  8kPa 
reduces  the  wave  speed  by  a  factor  of  ten.  This  also  reduces  Q  by  a  factor  of  ten. 18  Thus,  to  perform 
a  valid  experiment,  the  parameters  that  determine  Q  need  to  be  adjusted  in  the  softer  material  case  to 
compensate  for  the  slower  wave  speed  and  thus  restore  Q  to  a  larger  value. 

An  obvious  change  to  make  is  reducing  the  specimen  length.  Although  in  practice  this  might  be  quite 
difficult  for  brain  tissue,  reducing  the  length  (i.e.,  thickness)  of  the  specimen  by  a  factor  of  10  (from  1.7 

lsThe  symbol  Q  without  an  r  or  p  subscript  refers  to  both  Qr  and  Qp. 


_  lP 


Ctr, 


^travel 


(3-20) 
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to  0.17  mm)  would  increase  Q  by  a  factor  of  10.  However,  by  reducing  the  length,  inhomogeneities  in 
the  tissue  are  more  likely  to  influence  measurements.  There  are  also  practical  concerns:  how  to  cut  slices 
that  arc  submillimeter,  how  to  handle  them,  and  the  effects  of  gluing  them. 

Another  means  of  increasing  the  quality  Qr  is  to  increase  the  rise  time  of  the  loading  pulse.  To  increase 
Qr  by  a  factor  of  10  without  reducing  the  specimen  length,  tr  must  also  increase  by  a  factor  of  10;  i.e.,  we 
need  tr  =  2500  /is.  From  (3-19)i,  we  see  that  this  can  be  achieved  by  reducing  the  maximum  (average) 
strain-rate  ymax  or  by  increasing  the  maximum  (average)  strain  ymax.  However,  in  an  actual  torsion  test, 
buckling  or  tearing  is  likely  to  occur  at  sufficiently  large  strain.  Similar  comments  apply  to  Qp  and  tp. 

There  are  experimental  constraints  on  the  duration  of  the  loading  pulse,  namely  the  bar  length,  bar 
wave  speed,  and  location  of  the  strain  gage  on  the  bar.  Attempting  to  accommodate  a  long  (e.g.,  7  ms) 
loading  pulse  for  an  aluminum  or  steel  bar  by  only  increasing  its  length  would  yield  impractical  bar 
lengths.  However,  a  considerable  increase  in  loading  pulse  duration  can  be  gained  by  switching  to 
polymer  bar's  without  changing  the  length  of  the  bar-.  For  example,  polycarbonate  (PC)  and  polymethyl 
methacrylate  (PMMA)  have  shear-  wave  speeds  of  0.91  mm//zs  and  1.3  mm/gs,  respectively.  Compared 
to  aluminum  or  steel  bar's,  which  have  shear-  wave  speeds  on  the  order  of  3  nmi//xs,  the  longest  pulse  for 
a  fixed  length  of  bar-  would  be  increased  by  3.6  and  2.5  for  a  PC  and  a  PMMA  bar-,  respectively. 

Additional  gains  can  be  made  in  accommodating  longer  loading  pulses  by  replacing  the  strain  gage 
with  an  optical  measurement  of  the  angle  of  twist  imposed  on  the  specimen.  The  traditional  analysis  of 
Kolsky  bar  test  data  uses  the  fact  that  the  average  strain-rate  in  the  specimen  is  proportional  to  the  strain 
generated  in  the  incident  bar  by  the  reflected  wave,  assuming  the  specimen  is  in  dynamic  equilibrium ; 
see  (7)  in  [Gilat  2000]  for  torsion  tests  and  [Gray  2000;  Chen  and  Song  2011]  for  compression  tests.  For 
sufficiently  long  loading  pulses,  however,  the  incident  and  reflected  waves  can  overlap  at  the  strain  gage 
location,  in  which  case  the  reflected  wave  signal  cannot  be  separated  from  the  combined  signal,  making 
it  impossible  to  extract  the  strain  in  the  incident  bar-  associated  with  the  reflected  wave.19  This  would  no 
longer  be  an  issue  with  an  optical  measurement  of  the  angle  of  twist  at  the  end  of  the  incident  bar  adjacent 
to  the  specimen.  Also,  since  the  specimens  could  not  have  been  in  dynamic  equilibrium  in  the  tests  of 
[Nie  et  al.  2013],  the  method  discussed  above  for  determining  the  strain-rate  in  the  specimen  can  be  called 
into  question.  This  issue  would  also  be  bypassed  by  direct  optical  measurement  of  the  angle  of  twist. 

3F.  A  parametric  study  of  test  quality.  In  view  of  the  discussion  in  the  previous  section,  we  conducted 
a  parametric  study  in  which  tr  and  tp  were  varied  in  an  attempt  to  increase  the  quality  of  the  test  within 
the  experimental  constraints.  For  these  simulations,  we  used  ymax  =  700/s,  n  =  8kPa,  and  L  —  1.7  mm. 
With  these  parameters  fixed,  and  varying  only  the  rise  time  and  plateau  time,  we  present  in  Figure  5  the 
results  on  the  peak  strain  reached  at  the  Load  Cell  during  the  simulation  (panel  a),  the  total  duration  of 
the  loading  pulse  (panel  b),  and  the  average  error  (panel  c)  as  defined  below. 

Equation  (3-19)  predicts  that,  for  a  fixed  maximum  average  strain-rate  ymax ,  the  maximum  average 
strain  ymax  depends  linearly  on  the  sum  of  tr  and  tp.  In  panel  a ,  the  actual  peak  strain  in  the  specimen  at 
the  Load  Cell ,  max[y],  i.e.,  the  maximum  value  of  y(L,  t)  over  the  duration  of  the  experiment,  deviates 

19  This  was  one  of  the  difficulties  encountered  by  Nie  et  al.  [2013],  They  partially  bypassed  the  problem  by  substituting  the 
(negative  of  the)  incident  wave  for  the  reflected  wave,  arguing  that  this  is  a  good  approximation  for  very  soft  specimens,  which 
in  fact  it  may  be.  They  could  then  use  the  measured  strain  in  the  incident  bar  due  to  the  incident  wave  at  least  up  until  the  time 
of  arrival  of  the  reflected  wave. 
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Figure  5.  Panels  a-c  show  the  effects  of  varying  the  rise  time  tr  (horizontal  axis)  and  the 
plateau  time  tp  (vertical  axis)  for  ymax  =  700/s  tests  on  a  linear  elastic  specimen  with 
shear  modulus  8  kPa.  Through  false-color  images,  these  panels  show  the  peak  strain 
reached  at  the  Load  Cell  during  a  simulation  (max[y  ],  panel  a),  the  total  loading  pulse 
duration  tci  (panel  b),  and  the  average  error  (e)  in  the  strain  (see  text,  panel  c).  Panel  d 
shows  the  regions  where  the  peak  strain  is  less  than  1.5,  the  total  pulse  duration  less 
than  5  ms,  and  the  average  error  less  than  40%,  20%,  and  15%  in  darkening  shades 
of  gray.  The  black  dot  represents  Simulation  III,  which  most  closely  resembles  the  test 
conditions  in  [Nie  et  al.  2013].  The  red,  blue,  and  green  circles  are  sample  test  conditions 
taken  in  each  region.  Additional  results  for  these  test  conditions  are  given  in  panels  e 
and/.  Panel  e  is  a  plot  of  e(t)  (the  relative  error  in  the  strain  at  the  Load  Cell )  versus 
the  average  imposed  strain  ys(t).  Panel/  shows  the  resultant  stress-strain  curves.  The 
black  dashed  line  represents  the  ideal  to  which  these  curves  should  be  compared,  that  is, 
the  curve  that  would  have  been  produced  if  ys(t)  coincided  with  yLoadCeii- 


from  this  simple  behavior.  This  is  most  clearly  seen  for  the  smaller  values  of  tr  and  tp.  For  example, 
when  tr  =  250  /xs  and  tp  =  179  /xs,  the  maximum  average  strain  ymax  =  0.3,  but  the  peak  strain  is  close 
to  1  ( Simulation  III).  These  larger  strains,  as  discussed  in  the  previous  section,  are  from  a  failure  to 
achieve  dynamic  equilibrium. 

Panel  b  shows  the  dependence  of  the  pulse  duration  tc/  on  the  rise  time  and  plateau  time,  also  given 
by  (3-19).  This  panel  is  relevant  to  the  constraints  introduced  from  the  lengths  of  the  bars  in  a  Kolsky  bar 
experiment.  As  discussed  earlier,  a  7  ms  long  pulse  may  be  unmanageable  in  an  aluminum  or  steel  bar. 
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Since  there  are  nonuniform  strains  experienced  in  the  specimen,  there  is  an  error  associated  with 
assuming  that  the  measured  average  shear  strain  ys(t)  is  equal  to  the  shear  strain  at  the  Load  Cell,  y(L,  t). 
We  therefore  take  the  relative  error  in  the  strain  measurement  at  time  t  to  be 


_  I V(L,t)  -  ys(?) | 

Ks(0 


(3-21) 


This  eixor  is  time-dependent  and  is  small  when  near-uniform  conditions  exist  within  the  specimen.  Since 
y(L,t)  =  0  prior  to  the  arrival  of  the  shear  wave  at  the  Load  Cell,  the  value  of  the  error  e(t)  in  (3-21) 
is  1  for  all  times  less  than  the  travel  time.  We  ignore  the  error  during  these  early  times  since  we  hold 
the  travel  time  fixed  in  this  section.  We  define  an  average  error  by  using  the  time  average  of  its  value 
starting  after  the  travel  time  until  a  later  time  /max 


|  r  ^max 

{e)  = -  /  e(t)dt.  (3-22) 

hnax  firavel  J  ftravei 

Qualitatively  similar  results  were  found  for  various  values  of  fmax,  and  in  what  follows,  we  chose  a 
large  integration  window  where  ?max  =  15  ms  (roughly  twice  the  value  of  the  longest  pulse  duration 
considered).  Thus,  the  average  error  is  an  estimate  of  the  sensitivity  of  a  torsion  test  to  the  parameters 
of  the  loading  pulse. 

Panel  c  in  Figure  5  presents  the  dependence  of  the  average  error  (e)  on  the  rise  time  and  plateau  time. 
This  panel  clearly  shows  that  error  is  reduced  by  increasing  both  the  rise  time  and  the  plateau  time  as 
suggested  by  (3-20)  although  it  depends  most  critically  on  the  rise  time.  For  sufficiently  large  rise  times, 
the  effects  of  the  plateau  time  become  negligible. 

Combining  the  limitations  introduced  by  the  specimen  in  terms  of  maximum  allowable  strain,  the  con¬ 
straint  from  bar  lengths  in  terms  of  total  duration  of  the  pulse,  and  the  acceptable  error  in  the  measurement, 
we  can  determine  whether  it  is  possible  to  design  an  experiment  that  reaches  dynamic  equilibrium  for  a 
given  imposed  average  strain-rate.  The  small  dark  gray  region  in  panel  d  in  Figure  5  (small  triangular 
region  near  tr  =  2000  /xs)  represents  the  rise  and  plateau  times  where  the  peak  strain  at  the  Load  Cell  is 
less  than  1.5,  the  loading  pulse  duration  is  less  than  5  ms,  and  the  average  error  is  less  than  15%.  The 
two  neighboring  lighter  shades  of  gray  correspond  to  the  regions  where  rise  and  plateau  times  satisfy  the 
same  requirements  in  peak  strain  and  loading  pulse  duration  but  relax  the  average  error  constraint  to  20 
and  40%.  Although  not  clear  in  the  presentation  of  the  colored  regions,  these  regions  overlap  so  that 
the  light  gray  region  extends  “underneath”  the  two  darker  shades  of  gray,  and  the  medium  gray  region 
extends  under  the  dark  gray  region  (but  not  under  the  light  gray).  An  average  strain  of  1 .5  corresponds  to 
roughly  a  17°  rotation  of  one  end  of  the  specimen  with  dimensions  given  in  Table  2.  Achieving  a  pulse 
that  is  5  ms  long  would  require  a  polymer  bar  approximately  5  m  long,  depending  on  how  the  strain  was 
measured.  The  error  bounds  were  taken  fairly  loosely  since  actual  error  is  time-dependent  and  therefore 
depends  on  the  average  imposed  strain  ys.  The  red,  green,  and  blue  circles  in  panel  d  are  three  sample 
test  conditions,  one  taken  in  each  region  of  error.  The  black  dot  outside  the  gray  regions  represents 
Simulation  III,  which  most  closely  resembles  the  test  conditions  in  [Nie  et  al.  2013]. 

Panel  e  in  Figure  5  is  the  eixor  e(t)  plotted  against  the  imposed  average  strain  ys(t)  for  the  three  sample 
test  conditions  discussed  above.  This  panel  shows  that,  even  for  the  most  accurate  test  (green  line),  the 
initial  ring-up  period  extends  until  moderate  strains  of  about  0.4  are  reached.  This  means  that  even  in 
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these  experiments  the  small  and  moderate  strain  region  of  the  stress-strain  curve  would  be  unreliable. 
Panel  /  shows  the  resultant  stress-strain  curves  from  each  of  the  color-coded  simulations.  Here,  as  in 
Figure  4,  the  black  dashed  line  represents  the  ideal  to  which  these  curves  should  be  compared,  that  is,  the 
curve  that  would  have  been  produced  if  ys(t)  coincided  with  yu>ad  Cell ■  Long  loading  pulses20  and  large 
strains  are  necessary  to  reproduce  the  linear  stress-strain  curve  in  the  most  accurate  test  (green  line). 

4.  Quasistatic  torsion  of  isotropic  nonlinear  elastic  specimens21 

In  this  section,  we  examine  the  stress  state  in  isotropic  nonlinear  elastic  materials  undergoing  various 
types  of  quasistatic  torsional  deformation.  The  general  constitutive  relation  for  incompressible  isotropic 
elastic  solids  is  discussed  in  Section  4 A.  The  strain  state  in  pure  torsion  is  considered  in  Section  4B.  The 
stress  state  for  pure  torsion  of  a  general  incompressible  isotropic  elastic  solid  is  derived  in  Section  4C, 
and  the  stress  state  for  the  special  case  of  a  Mooney-Rivlin  material  is  considered  in  Section  4D.  The 
specimens  in  Sections  4B-4D  may  be  either  annular  or  disc-shaped.  For  incompressible  annular  speci¬ 
mens,  radially  nonuniform  torsional  deformations  are  briefly  discussed  in  Section  4E;  such  deformations 
arc  necessary  in  order  to  satisfy  the  stress-free  boundary  condition  on  both  the  outer  and  inner  surfaces. 

In  Section  4F,  we  present  the  results  of  our  three-dimensional  quasistatic  torsion  simulations  for  short, 
nearly  incompressible,  annular  specimens  with  glued  ends  and  stress-free  inner  and  outer  surfaces.  The 
constitutive  model  used  in  these  simulations  is  a  compressible  version  of  the  Mooney-Rivlin  model.  The 
boundary  conditions  in  these  simulations,  which  mimic  those  in  the  tests  by  Nie  et  al.  [2013],  do  not 
allow  the  annular  specimens  to  undergo  pure  torsion  or  the  radially  nonuniform  torsion  considered  in 
Section  4E.  The  stress  components  obtained  in  the  simulations  are  compared  with  the  analytical  results 
for  pure  torsion  from  Sections  4C  and  4D. 

4A.  Incompressible  isotropic  nonlinear  elastic  solids.  The  Cauchy  stress  tensor  T  can  be  decomposed 
into  a  pressure  (or  hydrostatic  stress)  p  and  a  deviatoric  stress  tensor  S  (i.e.,  tr  S  =  0): 

T  —  —pi  +  S,  p  =  -^trT,  S  =  devT  =  T  -  ^(trT)I.  (4-1) 

Here  I  denotes  the  identity  tensor.  The  standard  sign  convention  for  the  Cauchy  stress  tensor  is  used  in 
this  paper:  normal  stress  components  are  taken  positive  in  tension.  However,  the  pressure  p  is  positive 
in  compression.  It  is  sometimes  simpler  to  decompose  T  as 

T  =  -pi  +  5,  (4-2) 

A  A  r\r\ 

where  the  tensor  5  need  not  be  deviatoric,  in  which  case  the  scalar  p  need  not  equal  the  pressure  p. 
Consistency  with  the  relations  (4-1)  requires  that 

p  —  p—  |trS,  5  — dev  5.  (4-3) 


-°In  this  case,  the  notion  of  “constant  strain-rate”  should  be  modified  to  “peak  strain-rate”  since  the  plateau  time  is  small 
compared  to  the  rise  time. 

2 'General  background  for  some  of  the  material  in  Sections  4A-4E  can  be  found  in  the  books  [Green  and  Adkins  1960; 
Truesdell  and  Noll  1965;  Treloar  1975;  Batra  2006]  and  the  papers  [Rivlin  1948b;  1949;  Beatty  1987]. 

22  Nevertheless,  some  authors  still  refer  to  p  as  “the  pressure”  in  this  case. 
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For  a  compressible  (i.e.,  unconstrained)  material,  a  constitutive  relation  must  be  provided  for  T  and, 
equivalently,  for  p  and  5  or  for  p  and  S. 

For  an  incompressible  material,  p  and  p  are  indeterminate  in  the  sense  that  they  are  not  subject  to  any 
constitutive  relation.  However,  their  values  at  each  point  in  the  material  must  be  such  that  the  boundary 
conditions  and  momentum  balance  equations  (or,  in  the  quasistatic  case,  the  equilibrium  equations)  arc 
satisfied.  For  an  incompressible  material,  a  constitutive  relation  need  only  be  provided  for  5  or  5,  which 
arc  referred  to  as  the  determinate  parts  of  the  stress. 

For  an  incompressible  elastic  material,  the  determinate  stress  S  is,  by  definition,  a  function  of  the 
deformation  gradient  F .  If  the  material  is  also  isotropic,  then  S  is  an  isotropic  function  of  the  left 
Cauchy-Green  deformation  tensor  B  =  F Fr .  A  general  representation  for  such  functions  is  given  by 

5  =  n[(oB  -  (1  -  co)B~l]  =2C{B-  2 C2B~X ,  (4-4) 

where 

ii  Ci 

Ci  =  7jiuo,  C2  =  t/t(1  —  co),  p,  =  2(Ci  +  C2),  co  =  c  +  c  .  (4-5) 

The  elastic  moduli  Ci,  C2,  and  //  and  the  dimensionless  parameter  co  arc  scalar-valued  functions  of  tr  B 
and  tr  B~l .  In  particular,  /z  =  /z(tr  B  Ac  B  1 )  can  be  interpreted  as  a  (generally  strain-dependent)  shear 
modulus.23  For  infinitesimal  deformations,  the  deviatoric  paid  of  the  constitutive  relation  (4-4)  reduces 
to  the  linear  elastic  constitutive  relation  for  the  deviatoric  stress  tensor:  5  =  2/xq  dev  € ,  where  e  is  the 
infinitesimal  strain  tensor  and  /zo  =  /z( 3,  3)  since  tr  B  —  tr  B  1  =  tr  /  =  3  in  the  undeformed  state.  Thus, 
/zo  is  the  lineal-  elastic  shear  modulus  that  was  denoted  by  /z  in  previous  sections. 

Specification  of  either  set  of  constitutive  functions  /z,  co  or  Ci,  C2  determines  a  particular  material 
within  the  general  class  of  incompressible  isotropic  elastic  materials.  At  this  point,  we  do  not  assume 
that  the  material  is  hyperelastic  since  this  condition  is  not  necessary  for  any  of  the  results  on  the  stress 
state  in  pure  torsion  in  Section  4C.24  However,  the  incompressible  Mooney-Rivlin  model  considered 
in  Section  4D  and  the  compressible  version  of  this  model  used  in  the  three-dimensional  simulations 
(Sections  4F  and  5A)  are  hyperelastic. 

It  has  been  found  that  the  inequalities  C\  >  0  and  C2  >  0  imply  physically  reasonable  response  and 
are  consistent  with  experimental  data  on  many  nearly  incompressible  isotropic  elastic  solids.25  These 
inequalities,  which  are  equivalent  to 


/z>0,  0<<y<l,  (4-6) 

are  assumed  throughout  this  paper.  Then  0  <  1  —  co  <  1  also. 


23  See  also  (4-16)  below  for  the  case  of  pure  torsional  deformations. 

24  If  the  material  is  hyperelastic,  then  the  stress  tensor  is  derivable  from  a  strain  energy  (per  unit  volume)  W  with  W  = 
W(tr B,  ti -B~l)  for  an  incompressible  isotropic  elastic  material.  Consequently,  the  coefficient  functions  p,  co  and  C\,  C2  in 
(4-4)  are  determined  by  W',  e.g.,  C\  =  dW/d(tr B)  and  C2  =  3lV/3(trB_1).  It  follows  that  Cj  and  C2  are  not  independent 
although  neither  function  completely  determines  the  other.  Similar  comments  apply  to  /z  and  co  although  their  expressions  in 
terms  of  the  derivatives  of  W  are  more  complicated. 

25  See  [Truesdell  and  Noll  1965,  §49,  §53]  and  [Beatty  1987],  where  these  are  referred  to  as  the  E-inequalities. 
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4B.  Pure  torsional  deformation.  A  pure  torsional  deformation  satisfying  the  boundary  conditions  (2-2) 
is  given  by  (3-1)  with  the  angular  displacement  S  as  in  (3-13): 

r  —  R,  9  =  0(0,  Z)  =  0  +  t/r(Z-L),  ^  = z-  (4-7) 


We  assume  that  the  angle  of  twist  per  unit  length,  x/r,  is  positive:  i//  >  0.  For  this  deformation,  the 
physical  components26  of  the  deformation  gradient  F ,  the  left  Cauchy-Green  tensor  B  =  FFT ,  and  its 
inverse  B  1  arc  given  by 


o 

o 

o 

o 

o 

o 

0  1  t J/r 

[B]  = 

0  1  +  xj/2r2  xfrr 

[ B~l]  = 

0  1  — t J/r 

o 

o 

1 

V* 

O 

0  —  t jrr  1  +  i J/2r2 

(4-8) 


Thus,  pure  torsion  is  radially  inhomogeneous  and  isochoric  (det  F  =  1).  The  Eulerian  finite  strain  tensors 
corresponding  to  B  and  B  1  are  E  =  \(B  —  /)  and  £  =  \(I  —  B  ] );  the  latter  is  known  as  the 
Almansi  strain  tensor.  For  pure  torsion,  the  only  nonzero  components  of  E  and  £  are  Ege  =  f2r2, 
£zz  =  —%jr2r2,  and 

Eez  =  EZ0=Sez  =  Eze  =  \^rr,  (4-9) 


which  is  a  measure  of  the  shear  strain  in  pure  torsion.27 
The  volume  average  of  the  shear  strain  Egz  is  given  by 


(Eez) 


Eg-AV  =  \x!f 


R~  +  Rj  R0  +  R2 

Ri  +  Ro 


(4-10) 


Since  Egz  is  independent  of  z,  this  is  also  the  average  over  the  annular  cross-section  at  any  given  axial 
location.  (Egz)  can  also  be  expressed  in  terms  of  the  average  radius  R  and  the  radial  difference  A 
(see  (2-1)): 


(Eez)  =  & 


(3  R2  +  A2) 
2  R 


(4-11) 


where  the  approximation  is  valid  when  A/R  <£  1.  The  relative  error  in  approximating  Egz  by  ^xj/R  takes 
on  all  values  between  the  following  bounds: 

—  A/R  E0z-\xfR  A/R 

- - 2  (4-12) 

l  —  A/R  Eez  l  +  A/R 

The  relative  error  of  largest  absolute  value  occurs  at  the  inner  surface  and  is  given  by  the  absolute  value 
of  the  expression  on  the  left.  For  the  specimen  dimensions  considered  here  (see  Table  1),  the  largest 
relative  error  in  approximating  Egz  by  its  average  is  15%. 

As  discussed  in  Section  2B,  for  torsion  tests,  a  measure  of  average  shear  strain  often  used  in  the 
experimental  literature  is  ys  as  defined  in  (2-5),  that  is,  ys  =  f  R.  By  (4-11),  it  follows  that 


Ys~2{Egz)  (4-13) 

26  Recall  that  these  are  the  components  relative  to  the  unit  basis  vectors  along  the  cylindrical  coordinate  directions. 

27  In  general,  both  E  and  £  reduce  to  the  infinitesimal  strain  tensor  e  in  the  limit  of  small  displacement  gradients.  This  fact, 
together  with  (4-9),  yields  egz  =  ji/fr  for  infinitesimal  deformations,  consistent  with  (3-14). 
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with  the  relative  error  in  this  approximation  being  given  by 

2(Eez)  -  Ys  i  /  A\2 

Ks 


(4-14) 


For  the  specimen  dimensions  in  Table  1,  this  is  about  0.0056,  that  is,  a  relative  error  of  about  0.56%. 

Finally,  note  that,  while  annular  specimens  have  been  assumed  in  paid  of  the  discussion  above,  most 
of  the  results  in  this  section  are  also  valid  for  disc-shaped  specimens  ( R,  =  0).  The  approximation  on 
the  right  side  of  (4-11)  is  no  longer  valid  in  this  case  since  A/R  —  1.  Also,  the  lower  bound  in  (4-12) 
becomes  —  oo  since  A/R  ->  1  from  below  as  /?,  — >  0.  This  latter  result  also  follows  from  the  fact  that 
Egz  =  0  at  r  =  0  whereas  \  f  R  >  0. 


4C.  Stress  state  in  pure  torsion  of  incompressible  isotropic  elastic  specimens.  In  this  and  the  next 
section,  we  discuss  the  stress  state  in  pure  torsion  of  an  incompressible,  isotropic,  nonlinear  elastic  solid. 
The  results  in  these  sections  are  valid  for  both  annular  and  disc-shaped  specimens.  In  this  section,  we 
consider  the  general  constitutive  relation  (4-4)j  for  the  determinate  stress  S  subject  only  to  the  inequalities 
(4-6)  for  the  coefficient  functions  /x  and  a>.  In  the  next  section,  we  consider  the  special  case  where  these 
functions  are  taken  to  be  constants.  Recall  that  n  and  to  are  generally  scalar-valued  functions  of  tr  B 
and  tr  B  1 .  Since  tr  B  =  tr  B~l  =  3  +  x/rr2  for  pure  torsion,  we  have 

Ijl  —  /x(tr  B.  tr  B~l)  =  /x( 3  +  x/r2r 2 ,  3  +  i J/2r2)  =  jl(xj/2r2)  (4-15) 


in  this  case  and  similarly  for  a>,  Cj,  and  C 2- 

For  pure  torsion,  the  shear  stress  components  are  given  by  [Rivlin  1947;  1948b;  1949] 

Tdz  =  p.xj/r  =  2pEgz,  Trz  =  Trg  =  0,  (4-16) 

and  the  differences  in  the  normal  stress  components  are  given  by 

Tee  -  Tzz  =  p,\f2r2,  Tee  -  Trr  =  p,a)xf2r2,  Trr  -  Tzz  =  p(  1  -  co)\f2r2.  (4-17) 

These  relations  follow  from  (4-2),  (4-4)i,  and  (4-8).  It  is  clear  that  the  stress  state  is  radially  inhomoge¬ 
neous.  We  refer  to  the  shear  stress  component  77-  as  the  torsional  stress.  From  (4-16),  we  see  that  n 
can  be  interpreted  as  a  shear  modulus  and  that  Tgz  >  0  (since  we  assume  f  >  0)  with  the  exception  that 
Tgz  =  0  at  r  =  0  for  a  disc-shaped  specimen. 

From  the  relations  in  (4-16)  and  (4-17),  we  see  that  the  normal  stress  differences  arc  proportional  to 
the  torsional  stress: 


Tee  -  Tzz  =  fr  ■  Tgz,  Tee  -  Trr  =  cox/sr  ■  Tdz,  Trr  -Tzz-( l-  co)xjrr  ■  Tdz.  (4-18) 


The  relation  on  the  left  in  (4-18)  is  a  universal  relation  between  the  hoop,  axial,  and  torsional  stress:  it 
does  not  explicitly  involve  the  constitutive  functions  /i  or  o>  in  the  representation  (4-4)  i  for  the  determi¬ 
nate  stress. 

To  determine  the  normal  stress  components  (as  opposed  to  just  their  differences),  we  need  to  utilize 
the  equilibrium  equations  and  boundary  conditions.  Since  Trtj  and  Trz  arc  zero  throughout  the  specimen, 
the  radial  component  of  the  equilibrium  equation  reduces  to 


3  Trr  Trr  —  Tee  _  ^ 


dr 


r 


(4-19) 
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and  the  assumption  that  the  outer  lateral  surface  (. r  =  Ra)  is  stress- free  reduces  to  the  requirement  that 
Trr  =  0  at  r  =  R„.  On  using  (4-17)2  to  evaluate  the  Trr  —  Toe  term  above,  integrating  from  r  to  Ra,  and 
then  using  Trr \r=R0  =  0,  we  obtain  the  following  relation  for  the  radial  stress: 

r  R„  rRo 

Trr  =  —  f2  /  iuor  dr  =  —  2x//2  /  Cjrdr.  (4-20) 

Then  from  (4-17),  we  see  that  the  hoop  and  axial  stresses  are  given  by 

Toe  =  Trr  +  tuoxf2r2,  Tzz  =  Trr  -  /x(l  -  co)xl/2r2  (4-21) 

with  Trr  as  in  (4-20).  Since  /x  and  oo  are  functions  of  \f2r2,  it  follows  that  Trr,  Too,  Tzz,  and  Tgz  depend 
only  on  \f  and  r  while  Trtj  =  Trz  =  0.  These  results  imply  that  all  terms  in  the  6  and  z  components  of  the 
equilibrium  equation  are  identically  zero.  Thus,  pure  torsion  is  possible  in  any  incompressible,  isotropic, 
nonlinear  elastic  specimen  (disc-shaped  or  annular),  and  the  outer  surface  R  —  R()  may  be  taken  to  be 
stress-free  [Rivlin  1947;  1948b;  1949]. 28 

Since  /x,  o>  and  C\  are  functions  of  i J/2r2,  the  integral  for  Trr  cannot  be  evaluated  unless  C\  (equiv¬ 
alently,  the  product  /x&>)  is  specified.  However,  since  /x  and  co  are  positive  and  co  <  1  and  7/  ^  0,29 
important  qualitative  properties  of  the  normal  stress  components  may  be  inferred  directly  from  (4-20) 
and  (4-21),  as  summarized  below.30 

For  pure  torsion  of  an  incompressible,  isotropic,  nonlinear  elastic  specimen  (disc-shaped  or  annular) 
with  stress-free  outer  surface  r  =  Ra,  we  have: 

Conclusion  1.  Trr  is  negative  for  Rj  <  r  <  R().  Thus,  the  radial  stress  is  compressive  except  at  the  outer 
surface,  where  it  is  zero.  Furthermore,  the  magnitude  of  the  radial  stress  increases  with  increasing  radial 
distance  from  the  outer  surface. 

Conclusion  2.  Since  Trr  <  0  at  r  =  Rj,  for  annular  specimens  (where  /?,  >  0),  we  see  that  a  (nonzero) 
compressive  normal  stress  must  be  applied  on  the  inner  surface  in  order  that  radial  equilibrium  be  satisfied. 
If  this  compressive  stress  is  not  applied,  then  the  specimen  cannot  undergo  a  pure  torsional  deformation. 

Conclusion  3.  Since  Trr  <  0  for  R,  <  r  <  Ra,  it  follows  that  Tzz  <  0  for  Rj  <  r  <  Ra  and  Tzz  <  0 
at  r  =  R0  with  equality  holding  only  if  o>  =  I  at  r  =  Ra.  Thus,  the  axial  stress  is  compressive  everywhere 
except  possibly  at  the  outer  surface  r  =  R,,,  where  it  is  either  compressive  or  zero. 

Conclusion  4.  Since  Trr  =  0  at  r  =  R„,  it  follows  that  Too  \,  =r„  =  /xcm/f 2r2  >  0.  Thus,  the  hoop  stress  Too 
is  positive  (i.e.,  tensile)  at  the  outer  surface,  and  by  continuity,  Tee  must  be  tensile  for  some  radial  distance 
into  the  interior  of  the  specimen. 

Conclusion  5.  The  radial,  hoop,  and  axial  stresses  depend  on  both  constitutive  functions  /x  =  /x(i//2/'2) 
and  cu  =  6)(ij/2r2)  while  the  torsional  stress  is  independent  of  co.  It  follows  that  no  information  on  the  di¬ 
mensionless  constitutive  function  co  can  be  inferred  from  measurements  of  the  torsional  stress  Tez  alone.31 

“S  Rivlin  assumed  that  the  material  is  hyperelastic  although  that  assumption  is  not  necessary  for  this  conclusion  [Truesdell 
and  Noll  1965,  §57], 

Actually,  we  have  assumed  that  \//  is  positive,  but  the  sign  of  \jf  does  not  affect  Conclusions  1-5. 

30  Many  (but  not  all)  of  these  conclusions  can  be  found  in  the  references  cited  previously. 

3 1  This  conclusion  holds  even  if  we  assume  (as  is  reasonable)  that  the  material  is  hyperelastic  although  we  omit  the  proof  of 
this  result. 
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4D.  Pure  torsion  of  Mooney-Rivlin  materials.  The  classical  Mooney-Rivlin  model32  is  a  special  case 
of  the  general  constitutive  relation  for  an  incompressible  isotropic  elastic  solid.  It  is  obtained  by  taking 
the  coefficient  functions  C\  and  C2  (equivalently,  /x  and  (»)  in  the  expression  (4-4)  for  the  determinate 
stress  S  to  be  constants.33  The  special  case  C 2  =  0  (i.e.,  co  =  1)  is  known  as  a  neo-Hookean  material. 
For  infinitesimal  deformations,  the  deviatoric  part  of  the  Mooney-Rivlin  constitutive  relation  reduces 
to  5  =  2 /x  dev  e ,  where  e  is  the  infinitesimal  strain  tensor.  Thus,  the  constant  /x  in  the  Mooney-Rivlin 
model  agrees  with  the  shear  modulus  in  the  linear  elastic  theory  [Mooney  1940;  Rivlin  1947]. 

The  Mooney-Rivlin  model  is  one  of  the  simplest  isotropic  nonlinear-  elastic  constitutive  models,  and 
yet  it  still  exhibits  features  that  are  in  qualitative  agreement  with  experimental  data  on  many  soft,  nearly 
incompressible  materials  under  moderate  deformations.  For  this  model,  there  are  explicit  closed-form 
expressions  for  all  stress  components  in  pure  torsion  (see  below).  Furthermore,  a  compressible  version 
of  the  Mooney-Rivlin  model  was  available  in  the  code  that  we  intended  to  use  for  our  three-dimensional 
torsion  simulations.  For  these  reasons,  all  of  our  quasistatic  torsion  simulations  (Section  4F)  and  most  of 
the  three-dimensional  dynamic  simulations  (Section  5A)  used  this  compressible  version  of  the  Mooney- 
Rivlin  model.  The  following  analytical  results  for  pure  torsion  of  an  incompressible  Mooney-Rivlin 
material  will  be  compared  with  the  three-dimensional  quasistatic  torsion  simulations  in  Section  4F. 

For  Mooney-Rivlin  specimens,  the  relation  (4-16)  for  77-  and  the  relations  (4-17)  and  (4-18)  for  the 
normal  stress  differences  do  not  simplify  further  (aside  from  the  fact  that  /x  and  <»  are  now  constants). 
However,  the  relations  (4-20)  and  (4-21)  for  the  normal  stress  components  reduce  to34 

-Trr  =  \^2  ■  0){R2o  -  r2)  =  C\f2{R2  —  r2),  (4-22a) 

Tee  =  ■  <w(3r2  -  R2)  =  C^2(3r2  -  R2),  (4-22b) 

-Tzz  =  i/x^2  •  [coRl  +  (2  -  3 co)r2]  =  f2[{C,  -  2 C2)(R2  -  r2)  +  2C2R2].  (4-22c) 

Of  course,  all  of  the  conclusions  at  the  end  of  Section  4C  apply  to  the  special  case  of  Mooney-Rivlin 
specimens.  The  (compressive)  radial  stress  that  must  be  applied  to  the  inner  surface  r  =  7?,-  of  an  annular 
specimen  in  order  to  maintain  a  pure  torsional  deformation  is  now  given  explicitly  by 

Trr \r=R,  =  ■  0>{R2o  -  Rf)  =  ~f2Cx  (R2  -  R2)  <  0.  (4-23) 

From  (4-22b),  it  follows  that  Tee  >  0  (i.e.,  the  hoop  stress  is  tensile)  for  all  R,  <  r  <  R0  if  and  only 
if  Rj  >  R0/V 3.  This  condition  on  the  inner  and  outer  radii  is  satisfied  for  the  specimens  in  our  quasistatic 
and  dynamic  simulations  (see  Table  1).  The  relations  (4-22),  together  with  p  —  —  \{Trr  +  Tee  +  Tzz), 
imply  that  the  pressure  is  given  by 

p  =  \pLxlr2-[a)R20  +  \(2-l(o)r2].  (4-24) 

Since  /x  is  constant,  from  the  relation  (4-16)  for  Toz,  together  with  the  expression  (4-11)  for  {Eez)  and 


32  In  some  of  the  literature,  this  is  referred  to  as  the  Mooney  model  [Mooney  1940;  Rivlin  1947;  1948a;  1948b;  1949;  Rivlin 
and  Saunders  1951;  Green  and  Adkins  1960;  Truesdell  and  Noll  1965.  §95;  Treloar  1975;  Beatty  1987;  Batra  2006]. 

33  Equivalently,  this  is  the  hyperelastic  material  with  strain  energy  function  W  =  CjCti  B  —  3)  +  C2(trB_1  —  3). 

34  The  relations  (4-22c)2  and  (4-23)2  are  some  of  the  earliest  results  of  [Rivlin  1947;  1948a]. 
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the  fact  that  y^  —  tj/R  (see  (2-5)),  we  see  that  the  volume-averaged  torsional  stress  is  given  by 

(T0z)  =2p,{E9z)  =  p,fR  i  +  ld)2  1  +  Kf)2  •  (4'25) 

For  the  specimen  dimensions  in  Table  1,  the  term  in  square  brackets  differs  from  unity  by  only  0.0056. 
Thus,  the  relative  error  in  approximating  (7%-)  by  /xys  is  0.56%,  a  result  which  also  follows  directly 
from  (4-14). 

4E.  Radially  nonuniform  torsional  deformations  of  annular  specimens.  Conclusion  2  in  Section  4C 
implies  that  for  incompressible  annular  specimens  a  class  of  torsional  deformations  more  general  than 
pure  torsion  must  be  considered  if  both  the  outer  and  inner  lateral  surfaces  (r  =  R0  and  r  =  Rf)  are  stress- 
free.  Recall  that  the  pure  torsional  deformation  (4-7)  does  not  involve  any  radial  displacement.  Here 
we  briefly  consider  more  general  torsional  deformations  for  which  the  annular  specimen  may  deform 
radially: 

r  =  r(R ),  <9  =  0  +  VHZ-L),  z  =  Z.  (4-26a) 

Incompressibility  requires  that  1  =  dct  F  =  (dr/l) R)(r/ R)  =  3 (r2)/d(R2);  consequently, 

r  =  s/R2  +  P,  (4-26b) 

where  ft  is  an  arbitrary  constant  subject  only  to  the  condition  that  Rj  +  /3>  0.  Note  that  the  deformation 
(4-26)  reduces  to  pure  torsion  (r  =  R)  if  and  only  if  f  =  0. 

Proceeding  as  in  the  case  of  pure  torsion,  it  can  be  shown  that  the  deformation  (4-26)  is  possible  in  any 
incompressible,  isotropic  elastic,  annular  specimen;  that  Trz  =  Tro  =  0;  and  that  the  outer  surface  R  =  Ra 
may  be  taken  to  be  stress-free.  This  was  first  shown  by  Rivlin  [1949] . 35  The  radial  stress  is  expressed  in 
terms  of  the  undeformed  radial  coordinate  by 

r„  (S)  =  -  £  „[>«*  4-  i  -  d*  <4-27, 

with  the  constitutive  functions  //  and  o>  generally  dependent  on  i//2.  R2,  and  f.  On  letting  R  =  Rj  in  this 
relation  and  then  setting  the  result  to  zero,  we  obtain  an  implicit  relation  for  f>  that  must  be  satisfied  in 
order  that  the  inner  surface  R  =  R,  is  stress-free.  Clearly,  this  value  of  f>  depends  on  <//,  Rj.  and  R()  as 
well  as  the  constitutive  functions  /x  and  a>.  The  dependence  on  /x  drops  out  if  /x  is  a  constant  as  in  the 
Mooney-Rivlin  model.  For  that  case,  the  integral  can  be  evaluated  and  a  simpler  implicit  relation  for  ft 
is  obtained  [Rivlin  1949,  §8]. 36 

4F.  Three-dimensional  quasistatic  torsion  simulations.  Although  the  torsional  deformation  (4-26)  is 
of  some  theoretical  interest,  in  practice,  it  would  be  difficult  to  apply  a  torque  on  the  face  of  a  soft 
specimen  while  simultaneously  admitting  a  radial  deformation  on  that  face.  In  the  torsional  Kolsky  bar 
tests  in  [Nie  et  al.  2013],  the  faces  of  the  specimen  arc  glued  to  the  bar  and  the  load  cell,  and  analogous 

15  Rivlin  considered  a  slightly  more  general  deformation  that  includes  a  superposed  uniform  axial  extension.  He  assumed 
that  the  material  is  hyperelastic,  but  this  condition  is  not  necessary  [Truesdell  and  Noll  1965,  §57], 

36  Rivlin  considered  the  more  general  case  that  includes  a  superposed  uniform  axial  extension;  see  also  (95.8)  in  [Truesdell 
and  Noll  1965], 
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glued  end  conditions  are  used  for  their  quasistatic  torsion  tests.  Thus,  for  the  actual  torsion  tests  of 
interest  here,  neither  the  pure  torsion  considered  in  Sections  4B-4D  nor  the  torsional  deformation  (4-26) 
can  be  expected  to  hold  exactly.  Nevertheless,  it  is  clear  from  the  discussion  in  Section  4E  that,  at 
axial  locations  between  the  two  faces,  an  annular  specimen  will  have  a  tendency  to  deform  radially  in 
order  to  satisfy  the  stress-free  boundary  conditions  on  the  inner  and  outer  surfaces.  Since  relatively  short 
annular  specimens  are  used,  the  constraints  at  either  end  of  the  specimen  can  be  expected  to  severely 
restrict  (though  not  completely  eliminate)  the  radial  deformation  at  axial  locations  in  the  interior.  Thus, 
in  quasistatic  torsion  tests  with  glued  ends,  we  expect  some  axial  variation  in  the  radial  deformation  and 
consequently  in  the  stress  components  as  well  as  additional  nonzero  stress  components  such  as  Trz.  The 
analytical  results  for  the  torsional,  axial,  and  hoop  stresses  in  pure  torsion  may  still  hold  approximately, 
but  the  accuracy  of  these  approximations  is  not  clear.  Since  a  theoretical  analysis  of  this  problem  is 
intractable,37  we  rely  on  numerical  simulations  for  quantitative  estimates  of  the  stress  state. 

One  face  of  the  annular  specimen  was  subjected  to  a  rigid  rotation  while  the  other  face  was  fixed  so 
that  material  points  on  either  face  could  not  undergo  axial  or  radial  displacements.  Stress-free  boundary 
conditions  were  used  on  the  outer  surface  R  =  R„  and  the  inner  surface  R  =  Rj.  Details  of  the  numerical 
method  arc  provided  in  Appendix  A. 

In  these  simulations,  we  used  a  compressible  version  of  the  Mooney-Rivlin  model  for  the  specimens.38 
The  equation  of  state  for  the  pressure  in  this  model  is 


p  =  —  k  In  /  ~  —  J),  J  =  detF. 


(4-28) 


The  constant  k  is  the  bulk  modulus.  The  constitutive  relation  for  the  Cauchy  stress  tensor  T  is  given 
by  (4-2)  with  p  and  S  defined  as  follows: 

p  =  p+itr5,  JS  =  p[coB-(l-oo)B-l]  =  2ClB-2C2B~i.  (4-29) 


Here  /z,  co,  Cj,  and  C2  arc  constants  as  in  the  incompressible  version,  and 


B  = 


B 

(det  B)1/3 


B 

J2j3 


(4-30) 


is  the  distortional  part  of  the  left  Cauchy-Green  tensor  B.  Note  that,  regardless  of  the  value  of  the 
Jacobian  J,  det  B  =  det  B  1  =  1,  analogous  to  the  conditions  det  B  =  det  B~l  =  1  for  the  incompressible 
case.  This  relation  for  S  is  identical  to  the  relation  (4-4)  for  the  incompressible  Mooney-Rivlin  model 
except  that  B  has  been  replaced  by  its  distortional  part  and  S  is  multiplied  by  the  factor  J.39  The 
relation  (4-29)i  for  p  guarantees  consistency  with  (4-3)i . 

We  used  k  =  2.3  GPa  for  all  simulations  and  varied  the  shear  modulus  by  almost  three  orders  of  mag¬ 
nitude:  p,  =  800,  80,  8,  and  2kPa.  The  corresponding  ratios  of  bulk  to  shear  modulus  varied  from  2875 
to  1.15  x  106  so  that  in  all  cases  the  specimen  was  nearly  incompressible.  For  each  of  the  four  values 


33  Some  idea  of  the  difficulty  of  this  problem  can  be  inferred  from  the  paper  [Hill  and  Lee  1989],  which  treats  a  slightly 
different  problem,  namely  combined  compression  and  torsion  of  incompressible,  disc-shaped,  Mooney-Rivlin  specimens  with 
glued  ends.  This  work  has  been  extended  to  torsional  waves  [Hill  and  Wegner  2004],  but  the  results  are  not  directly  applicable 
to  the  dynamic  problem  considered  here  since  annular  specimens  are  not  considered. 

38  Additional  departures  (perhaps  small)  from  pure  torsion  will  be  introduced  by  the  slight  compressibility  of  the  specimens. 

39  This  factor  ensures  that  the  material  is  hyperelastic.  The  strain  energy  function  is  given  by  W  =  k  J (\n  J  —  \)  + 
Crftrfi  — 3)  +  C2(triH  -3). 
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of  the  shear  modulus,  simulations  were  performed  for  three  values  of  the  nondimensional  parameter  co: 
0.3,  0.6,  and  1.  This  resulted  in  twelve  quasistatic  simulations. 

We  have  compared  the  stress  states  obtained  from  the  numerical  simulations  (using  the  compressible 
version  of  the  Mooney-Rivlin  model)  with  the  stress  states  that  would  have  been  obtained  in  pure  torsion 
of  an  incompressible  Mooney-Rivlin  material  (as  determined  from  (4-22)  and  (4-24)).  These  compar¬ 
isons  are  made  at  the  same  average  angle  of  twist  per  unit  length  i/z  (equivalently,  at  the  same  average 
shear  strain  ys )  and  for  the  same  set  of  twelve  values  for  /z  and  co.  The  results  are  plotted  in  Figure  6. 
Each  of  the  ten  panels  contains  either  thirteen  or  fifteen  curves,  but  many  (in  panels  a  and  b,  all  of  them) 
are  indistinguishable.  This  is  due  (in  part)  to  the  normalization  of  the  stress  components  as  discussed  in 
the  next  paragraph.  In  each  panel,  there  are  twelve  curves  for  the  numerical  simulations  (solid  lines)  and 
either  one  curve  (panels  a  and  b)  or  three  curves  for  the  theoretical  pure  torsion  (dashed  lines). 


Figure  6.  Normalized  stress  states  in  torsion  for  Mooney-Rivlin  specimens.  Twelve 
numerical  simulations  (solid  lines)  are  compared  with  pure  torsion  (dashed  lines).  For 
each  of  four  values  of  the  shear  modulus  /z  (800,  80,  8,  and  2  kPa),  curves  for  co  =  0.3 
(red),  co  =  0.6  (blue),  and  co  =  1.0  (black)  are  shown.  Curves  for  different  shear  moduli 
are  indistinguishable  (either  approximately  or  exactly)  because  of  the  normalizations. 
The  left  panels  examine  the  variation  of  the  normalized  stress  with  referential  radius  for 
an  angle  of  twist  per  unit  length  i/r  =  0.059  radians/mm  at  the  axial  location  Z  =  L/2. 
The  right  panels  are  plots  of  normalized  volume-averaged  stress  versus  average  shear 
strain.  The  vertical  line  in  the  right  column  indicates  the  average  shear  strain  for  the 
panels  on  the  left. 
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In  order  to  condense  the  information  from  the  simulations  and  the  theoretical  results  and  also  to  more 
clearly  display  the  relative  magnitude  of  the  stresses  in  each  case,  all  stress  components  in  Figure  6  have 
been  normalized  by  an  appropriate  “maximum”  torsional  stress.  For  the  left  column,  the  tilde  on  the 
stresses  indicates  that  they  have  been  normalized  by  fix//  R0,  which  by  (4- 1 6)  i  is  the  maximum  torsional 
stress  that  would  occur  in  a  pure  torsional  deformation.  For  the  right  column,  the  hat  on  the  stresses 
indicates  that  they  have  been  normalized  by  the  volume-averaged  torsional  stress  (Tgz)  (see  (4-25))  that 
would  occur  in  a  pure  torsional  deformation,  evaluated  at  the  final  average  shear  strain  ys  =  0.5.  The 
angle  brackets  on  the  stresses  in  the  right  column  indicate  that  these  normalized  stresses  have  then  been 
volume-averaged. 

In  the  left  column  of  Figure  6,  we  examine  the  radial  variation  in  the  torsional  stress,  the  normal  stress 
components,  and  the  pressure  at  an  axial  location  in  the  middle  of  the  specimen  (Z  =  L/2).  In  these  plots, 
the  average  angle  of  twist  per  unit  length  is  x[r  =  0.059  radians/mm  =  3.38°/mm,  and  the  corresponding 
average  shear  strain  is  ys  =  0.38.  Data  from  torsion  tests  is  typically  plotted  as  average  shear  stress 
versus  average  shear  strain  ys  with  the  average  shear  stress  inferred  from  the  measured  torque.  The  right 
column  of  Figure  6  displays  the  volume-averaged  stress-strain  curves.  The  vertical  black  lines  in  the 
right  column  indicate  the  average  shear  strain  (0.38)  for  the  panels  on  the  left. 

Consider  panel  a  for  the  radial  variation  of  the  torsional  stress.  Since  Tgz  =  iix/jr  for  pure  torsion 
of  an  incompressible  Mooney-Rivlin  material,  the  normalized  torsional  stress  in  this  case  is  given  by 
Tez  =  fix/zr/ iixj/R0  =  r/Ra  =  R/R„,  independent  of  /x  or  co.  Thus,  for  pure  torsion  of  an  incompressible 
Mooney-Rivlin  material,  there  is  a  single  normalized  torsional  stress  versus  radius  curve,  namely  the 
line  R/R0.  This  line  and  the  twelve  normalized  torsional  stress  versus  radius  curves  from  the  numerical 
simulations  arc  indistinguishable,  the  maximum  relative  difference  being  less  than  1%. 

A  similar  situation  occurs  for  the  average  torsional  stress  versus  average  shear  strain  curves  in  panel  b. 
Since  (Tgz)  is  given  by  (4-25)  for  pure  torsion  of  an  incompressible  Mooney-Rivlin  material,  the  nor¬ 
malized  average  torsional  stress  in  this  case  is  given  by  (Tgz)  =  ys /(ys I o.s) -  Thus,  for  pure  torsion  of 
an  incompressible  Mooney-Rivlin  material,  there  is  a  single  normalized  average  torsional  stress  versus 
average  strain  curve,  namely  the  line  ys/0.5.  This  line  and  the  twelve  normalized  average  torsional  stress 
versus  average  shear  strain  curves  from  the  numerical  simulations  are  indistinguishable,  the  maximum 
relative  difference  again  being  less  than  1%. 

Thus,  even  though  the  annular  specimens  in  our  numerical  simulations  cannot  undergo  a  pure  torsional 
deformation,  the  torsional  stress  in  those  simulations  is  indistinguishable  from  the  torsional  stress  in 
theoretical  pure  torsion  of  an  incompressible  Mooney-Rivlin  material.  On  the  other  hand,  it  is  clear 
from  Figure  6  that  this  conclusion  does  not  hold  for  the  normal  stress  components;  that  is,  the  normal 
stresses  and  the  pressure  in  the  annular  specimen  simulations  (solid  lines)  do  not  agree  with  the  theoretical 
predictions  for  pure  torsion  of  an  incompressible  Mooney-Rivlin  material  (dashed  lines). 

For  pure  torsion  of  an  incompressible  Mooney-Rivlin  material,  the  theoretical  relations  (4-22)  and  (4-24) 
predict  that  the  normal  stress  components  and  the  pressure  are  proportional  to  the  shear  modulus  /x  and 
depend  on  co  as  well;  consequently,  the  normalized  stress  components  in  panels  c-j  in  Figure  6  are 
independent  of  /jl  but  still  depend  on  co.  This  results  in  three  normalized  stress  curves  (dashed  lines)  for 
the  case  of  pure  torsion  in  each  of  these  panels  —  one  curve  for  each  of  the  three  values  of  co. 

The  disagreement  between  the  normalized  stresses  for  pure  torsion  and  those  for  the  numerical  sim¬ 
ulations  are  most  pronounced  for  the  radial  component  (panels  c  and  cl).  In  fact,  the  normalized  radial 
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stresses  for  the  simulations  are  so  small  that  all  twelve  curves  (solid  lines)  are  indistinguishable  from 
the  horizontal  axis  (Trr  =  0  or  (Trr)  =  0).  Of  course,  the  fact  that  the  radial  stress  is  essentially  zero 
throughout  the  annular  specimen  is  not  unexpected  in  view  of  the  fact  that  it  is  zero  on  the  inner  and 
outer  lateral  surfaces.  This  is  one  of  the  motivations  for  using  annular  specimens.  The  other  motivation 
for  using  annular  (as  opposed  to  disc-shaped)  specimens  in  torsion  tests  is  that,  while  torsional  stress  Tez 
and  the  shear  strain  Egz  arc  radially  nonuniform,  the  total  radial  variation  in  these  quantities  is  expected 
to  be  smaller  for  annular  specimens.  For  the  torsional  stress,  this  is  verified  in  Figure  6 a. 

Unlike  the  radial  stresses,  the  hoop  stress,  axial  stress,  and  pressure  (panels  e—j)  for  the  simulations 
(solid  lines)  share  some  qualitative  features  with  the  theoretical  predictions  for  pure  torsion  (dashed  lines). 
In  particular,  the  normalized  hoop  stress,  axial  stress,  and  pressure  for  the  simulations  are  (approximately) 
independent  of  the  shear  modulus  but  depend  strongly  on  the  value  of  co  so  that  the  twelve  curves  appear 
to  collapse  into  three  curves.  Furthermore,  for  both  the  simulations  and  the  pure  torsion  case,  the  hoop 
stress  increases  with  increasing  co  while  the  magnitude  of  the  axial  stress  decreases  with  increasing  co. 
For  the  axial  and  hoop  stresses,  the  discrepancies  between  the  simulations  and  the  pure  torsion  case 
increase  with  increasing  co. 

Further  consideration  of  (4-16)  and  (4-22)  shows  that,  for  pure  torsion,  the  ratio  of  normal  stress  to 
torsional  shear  stress  increases  with  shear  strain:  the  shear  stress  depends  linearly  on  /  while  the  normal 
components  depend  quadratically  on  t,//  and  thus  quadratically  on  the  average  shear  strain  ys.  Comparison 
of  the  solid  curves  in  panel  b  with  those  in  panels/  and  /?  in  Figure  6  shows  that  the  same  conclusion 
holds  for  the  hoop  and  axial  stresses  (but  not  the  radial  stress)  in  the  numerical  simulations.  Thus,  for 
sufficiently  large  angles  of  twist  or  for  sufficiently  thin  specimens,  the  magnitude  of  the  normal  stresses 
could  exceed  the  shear  stress.  On  the  other  hand,  a  smaller  shear  strain  of  0.2,  say,  gives  an  average  hoop 
stress  that  is  only  5  to  20%  of  the  shear  stress  (depending  on  co)  and  an  average  axial  stress  that  is  only 
1  to  15%  of  the  shear  stress. 

The  glued-end  boundary  conditions  in  combination  with  the  stress-free  boundary  conditions  on  the 
inner  and  outer  surfaces  of  the  annular  specimen  resulted  in  axial  variations  in  the  stress  components  in 
the  torsion  simulations  although  these  results  are  not  plotted  here.  For  the  normal  stresses,  these  axial 
variations  were  concentrated  within  approximately  0.25  mm  of  either  face  of  the  specimen  and  essentially 
vanished  away  from  the  specimen  faces.  Both  hoop  and  axial  stresses  varied  by  as  much  as  30%  near 
the  faces  whereas  the  variation  in  this  region  for  the  torsional  stress  Tgz  was  typically  less  than  5%.  The 
other  shear  components,  Trz  and  Trfj.  were  also  nonzero.  They  were  smaller  in  magnitude  (typically  less 
than  10%  of  the  torsional  stress)  and  varied  throughout  the  length  of  the  specimen. 

Although  the  magnitude  of  Trz  is  small,  the  fact  that  Trz  cannot  be  identically  zero  can  be  established 
as  follows.  Using  the  assumption  that  the  specimen  is  in  equilibrium  (so  that  div  T  =  0)  together  with  the 
stress-free  boundary  conditions  on  the  inner  and  outer  surfaces  of  the  annular  specimen,  one  could  show 

f  -TggdA  =  ^-  f  TrzdA,  (4-31) 

J A(z)  r  oz  J _4(z) 

where  A(z)  is  the  annular  cross-section  at  the  axial  location  z.40  As  observed  above  (see  panel  e  in 


40  We  omit  the  proof  of  this  relation.  It  is  a  fairly  general  result  in  the  sense  that  no  material  symmetry  conditions  are 
required,  the  specimen  need  not  be  homogeneous,  and  neither  the  initial  nor  deformed  annular  cross-sections  need  be  axially 
uniform  or  have  circular  boundaries. 
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Figure  6),  Tee  is  positive  at  each  radial  location.  Hence,  the  integral  on  the  left  above  is  positive  and 
must  be  balanced  by  an  axial  gradient  in  the  integral  of  Trz  over  the  cross-section  A(z). 


5.  Three-dimensional  dynamic  simulations 

In  this  section,  we  further  supplement  our  one-dimensional  dynamic  simulations  and  three-dimensional 
quasistatic  analysis  and  simulations  with  three-dimensional  dynamic  simulations.  Our  dynamic  torsion 
simulations  use  the  same  glued-end  boundary  conditions  and  stress-free  boundary  conditions  on  the  outer 
and  inner  surfaces  as  in  the  quasistatic  simulations  discussed  in  Section  4F.  With  the  exception  of  the 
viscoelastic  results  discussed  in  Section  5B,  we  used  the  same  constitutive  model  and  the  same  model 
parameters  as  in  the  three-dimensional  quasistatic  simulations  except  that  the  /x  =  2  kPa  case  is  not  con¬ 
sidered.  Thus,  we  used  k  =  2.3  GPa  and  varied  the  shear  modulus  by  two  orders  of  magnitude:  n  =  800, 
80,  and  8  kPa.  The  corresponding  ratios  of  bulk  to  shear  modulus  varied  from  2875  to  2.875  x  105  so  that 
in  all  cases  the  specimen  was  nearly  incompressible.  For  each  of  the  three  values  of  the  shear  modulus, 
simulations  were  performed  for  three  values  of  the  nondimensional  parameter  co:  0.3,  0.6,  and  1.  The 
nine  simulations  used  the  same  700/s  loading  pulse  as  the  one-dimensional  Simulations  /-///;  see  Table  2 
and  panels  a  and  b  in  Figure  3.  Details  of  the  numerical  method  arc  provided  in  Appendix  A. 


5A.  Nonlinear  elastic  model  for  the  specimen.  We  calculated  the  radially  averaged  stress  at  the  Load 
Cell  for  each  of  the  six  components  T,}  of  the  stress  tensor,  i.e., 


1 


Tjj{R.  L)2izRdR. 


These  time  histories  arc  plotted  for  the  nine  simulations  in  Figure  7.  The  shear  modulus  /i  changes  with 
the  row,  and  co  changes  with  the  column  so  that  the  nine  panels  span  vastly  different  elastic  responses. 

In  the  quasistatic  torsion  simulations  (see  Figure  6),  the  torsional  shear  stress  Tqz  is  insensitive  to  the 
dimensionless  parameter  co  in  the  Mooney-Rivlin  model.  This  is  also  observed  in  the  dynamic  case: 
comparing  the  shear  stress  Tez  (green  traces)  for  fixed  shear  modulus,  i.e.,  panels  a-c,  d-f,  and  g-i 
in  Figure  7,  we  see  little  to  no  variation  with  co.  Similarly,  in  both  the  quasistatic  and  dynamic  simula¬ 
tions,  the  magnitude  of  the  axial  stress  Tzz  (black  traces)  decreases  with  increasing  co  whereas  the  hoop 
stress  Tee  (blue  traces)  increases  with  increasing  co  with  the  exception  of  the  hoop  stress  for  the  /r  =  8  kPa 
dynamic  case.  Since  neither  77,  or  Tee  is  measured  in  a  Kolsky  bar  test,  there  would  be  no  way  to  infer  a 
value  for  co  from  the  test  data,  namely  the  measured  shear  stress,  even  if  it  were  known  that  the  specimen 
behavior  was  governed  by  the  Mooney-Rivlin  model. 

Unlike  the  quasistatic  simulations,  for  a  given  value  of  co.  the  relative  magnitude  of  the  normal  and 
torsional  shear  stress  components  vary  considerably  with  the  shear  modulus  in  the  dynamic  simulations. 
In  particular,  the  radial  stress  (red  traces)  is  not  negligible  for  the  fi  =  80  and  8  kPa  dynamic  simulations. 
Also  note  that,  for  the  /i  =  8kPa  case  with  co  =  0.3  or  0.6  (panels  g  and  h),  all  three  normal  stress 
components  begin  to  change  well  before  the  torsional  stress  increases.  This  effect  cannot  be  due  to 
bending  waves  or  longitudinal  waves  that  might  be  generated  in  the  incident  bar  as  a  result  of  the  initial 
impact  since  such  effects  are  not  included  in  the  simulations.  Rather,  it  is  caused  by  a  longitudinal 
precursor  wave  in  the  specimen  as  discussed  in  the  next  paragraph.  From  Section  3,  it  is  clear  that 
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Figure  7.  The  radially  averaged  stress  components  at  the  Load  Cell  (Z  =  L)  versus 
time  for  the  nine  simulations:  two  decades  of  shear  modulus  (/x  =  800,  80,  and  8  kPa) 
and  three  values  of  dimensionless  parameter  (a>  =  0.3,  0.6,  and  1.0).  All  simulations 
used  the  same  700/s  loading  pulse  as  the  one-dimensional  Simulations  I -I II]  see  Table  2 
and  panels  a  and  b  in  Figure  3.  The  legend  for  the  line  colors  is  given  in  panel  /  and  is 
consistent  throughout  the  figure.  Note  the  change  in  vertical  scales  for  the  different  rows. 


dynamic  equilibrium  is  not  reached  when  /x  <  80 kPa,  and  the  three-dimensional  simulations  inform  us 
that  the  stress  state  is  no  longer  simple  as  well. 

The  normal  stresses  produced  by  the  shearing  deformation  propagate  at  the  longitudinal  elastic  wave 
speed,41  which  far  exceeds  the  shear  wave  speed  for  all  cases  in  Figure  7.  Indeed,  the  longitudinal  wave 
speed  is  about  1.52mm//xs  for  all  nine  cases.42  Flence,  this  wave  arrives  at  the  Load  Cell  in  about  1.1  /xs, 

41  This  is  a  well-known  nonlinear  elastic  effect  [Davison  1966;  Abou-Sayed  and  Clifton  1976;  Scheidler  1998], 

42  For  a  nonlinear  elastic  solid,  a  longitudinal  wave  propagating  into  undeformed  material  travels  at  the  linear  elastic  wave 
speed  cL  =  \l ^c|  +  c^,  where  cs  is  the  shear  wave  speed  and  cR  is  the  bulk  wave  speed,  i.e.,  *JKI P-  F°r  the  bulk  and  shear 
moduli  considered  here  (k  =  2.3  GPa  and  /x  <  800  kPa),  the  shear  wave  speed  is  negligible  and  cB  =  1.52  mm//xs. 
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which  is  essentially  instantaneous  on  the  time  scale  of  Figure  7.  But  even  though  the  normal  stresses 
begin  to  change  at  the  Load  Cell  at  this  early  time,  their  magnitude  is  initially  quite  small  —  on  the 
order  of  the  square  of  the  shear  strain  produced  by  the  slower  torsional  wave.  Consequently,  these  initial 
changes  cannot  be  observed  on  the  scale  of  the  plots  until  much  later  times.  And  it  is  only  in  panels  g 
and  h  that  the  magnitude  of  the  normal  stresses  (relative  to  the  torsional  stress)  is  significant  enough  to 
be  observed  prior  to  the  arrival  of  the  torsional  wave. 

5B,  Nonlinear  viscoelastic  model  for  the  specimen.  Although  the  goal  of  the  torsional  Kolsky  bar  test 
considered  here  is  to  quantify  the  strain-rate  sensitivity  of  viscoelastic  soft  tissues,  only  elastic  constitu¬ 
tive  models  have  been  considered  up  to  this  point.  Elastic  models  (linear  or  nonlinear)  are  particularly 
useful  for  revealing  inertial  effects:  since  the  stress  in  an  elastic  material  is  rate-independent,  any  dif¬ 
ferences  between  the  quasistatic  and  dynamic  simulations  of  the  torsion  test  must  necessarily  be  due 
to  inertial  effects.  There  is  certainly  no  reason  to  expect  that  these  inertial  effects  will  disappear  for 
materials  exhibiting  viscoelastic  response.  Nevertheless,  we  have  also  performed  some  three-dimensional 
dynamic  simulations  of  torsional  Kolsky  bar  tests  using  a  nonlinear  viscoelastic  constitutive  model43  for 
the  specimen. 

We  used  an  equilibrium  elastic  shear  modulus  of  2  kPa  and  an  instantaneous  elastic  shear  modulus 
of  40kPa,44  five  times  larger  than  the  8  kPa  shear  modulus  that  had  been  used  as  representative  of  brain 
tissue  in  the  elastic  models.  This  gives  a  viscoelastic  shear  wave  speed  about  2.2  times  the  shear  speed 
for  the  fi  =  8  kPa  linear  and  nonlinear  elastic  cases.45  It  should  be  clear  from  our  results  for  these 
elastic  cases  that  this  wave  speed  is  much  too  slow  for  dynamic  equilibrium  to  be  achieved.  Indeed, 
while  the  stress  histories  in  our  simulations  with  the  viscoelastic  model  differed  both  qualitatively  and 
quantitatively  from  those  described  above  for  the  elastic  model,  we  found  that  the  viscoelastic  stress 
histories  fell  roughly  between  the  n  =  80kPa  and  8  kPa  elastic  cases  in  Figure  7. 

6.  Discussion  and  conclusions 

We  have  conducted  one-  and  three-dimensional  analyses  and  simulations  to  understand  the  deformations 
and  stress  states  that  exist  in  nearly  incompressible  soft  materials  in  a  torsional  Kolsky  bar  test.  Here  we 
highlight  the  main  conclusions  of  this  study  and  also  discuss  some  relevant  issues  that  were  not  addressed 
previously. 

The  use  of  {linear  or  nonlinear)  elastic  constitutive  models  for  the  specimen  in  numerical  simulations 
of  Kolsky  bar  tests  can  reveal  the  presence  of  inertial  effects  in  these  tests.  Using  viscoelastic  models  in 

43  See  Appendix  A  for  more  details. 

44  This  is  close  to  our  rough  estimate  of  42kPa  for  the  700/s  tests  in  [Nie  et  al.  2013];  see  the  discussion  at  the  end  of 
Section  3D.  Recall  that  they  estimated  an  initial  shear  modulus  of  53.5  kPa  from  their  average  “700/s”  shear  stress  versus  shear 
strain  curve. 

45  For  a  (linear  or  nonlinear)  viscoelastic  material,  the  speed  of  a  shear  wave  propagating  into  an  undeformed  region  is  that 
of  a  linear  elastic  material  with  shear  modulus  equal  to  the  (small  strain)  instantaneous  elastic  shear  modulus  of  the  viscoelastic 
material.  This  result  also  holds  for  any  simple  material  with  fading  memory  [Coleman  et  al.  1965].  Flowever,  the  shear 
wave  in  the  torsion  test  initially  propagates  into  a  slightly  deformed  region  caused  by  the  longitudinal  precursor  wave  and 
then  reflects  from  the  Load  Cell  into  a  more  complicated  strain  state.  Nevertheless,  the  shear  wave  speed  into  an  undeformed 
material  should  yield  a  rough  estimate  for  the  wave  speed  in  this  test,  and  for  the  viscoelastic  model  considered  here,  this  is 
U40kPa/p  =  s/5  x  8  kPa / p  or  s/5  times  the  shear  speed  for  the  /x  =  8  kPa  elastic  cases. 
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these  simulations  will  add  strain-rate  effects  to  the  stress  histories  but  not  substantially  alter  the  inertial 
effects  since  the  main  issue  here  is  the  slow  shear  wave  speed  in  the  specimen. 

In  soft  specimens,  the  travel  time  of  the  torsional  wave  can  be  on  the  order  of  the  loading  pulse 
duration,  in  which  case  there  is  insufficient  time  for  ring-up  to  dynamic  equilibrium.  For  the  torsional 
Kolsky  bar  tests  on  bovine  brain  tissue  at  a  strain-rate  of  “700/s”  that  were  reported  in  [Nie  et  al.  2013], 
the  specimen  could  not  have  been  in  dynamic  equilibrium,  contrary  to  statements  made  in  that  paper. 
Consequently,  their  “stress-strain  curves”  are  not  valid  —  the  average  strain  used  in  their  plots  does  not 
correspond  to  the  strain  at  the  load  cell.  In  fact,  throughout  much  of  the  test,  the  actual  strain  at  the  load 
cell  was  likely  several  times  larger  than  the  average  strain  used  in  their  stress-strain  plots. 

The  stress  history  at  the  load  cell  should  never  be  time-shifted  to  synchronize  it  with  the  history  of 
the  average  strain  in  order  to  produce  a  “stress-strain  curve  ”.  If  the  stress-strain  curves  produced  from 
unshifted  data  contain  an  initial  strain  interval  for  which  the  stress  is  zero  (as  in  the  left  panels  in  Figure  4), 
then  this  is  a  clear  indication  that  the  specimen  is  not  in  dynamic  equilibrium  and  consequently  that  no 
meaningful  stress-strain  curves  can  be  inferred  from  the  test  data.  Even  if  the  “stress-strain  curves” 
produced  by  time-shifting  the  load  cell  data  appear  physically  reasonable,  they  will  necessarily  contain 
erroneous  “strain-rate  effects”. 

It  is  not  clear  how  to  modify  the  experimented  design  to  alleviate  the  inertial  effects  that  exist  when 
testing  with  very  soft  specimens  at  high  strain-rates.  In  Section  3E,  we  discussed  the  changes  in  the  load¬ 
ing  pulse  and  the  specimen  geometry  that  would  be  necessary  to  improve  the  quality.  These  alterations 
require  longer  bars,  thinner  specimens,  and  reductions  of  the  maximum  strain-rates  considered.  Ignoring 
the  practical  difficulties  that  longer  bars,  polymer  bars,  or  thinner  specimens  pose,  the  lower  strain-rates 
required  would  be  on  the  order  of  100-200/s  (unless  larger  strains  arc  considered)  and  are  perhaps  too 
low  for  the  ranges  needed  for  the  modeling  effort. 

Other  Kolsky  bar  tests  such  as  double-lap  shear  [Saraf  et  al.  2007;  Trexler  et  al.  2011]  or  circular 
shearing  [Nie  et  al.  201 1]  will  have  similar  issues  with  dynamic  equilibrium.  Similar  dimensional  argu¬ 
ments  can  be  made  for  each  of  these  tests.  With  proper  consideration  of  the  direction  of  the  shear  wave 
propagation  and  the  relevant  specimen  length,  the  results  of  Section  3  would  apply.  Thus,  inertial  effects 
are  still  an  issue  with  high  shear  strain-rate  Kolsky  bar  tests  on  soft  materials. 

In  view  of  the  remarks  above,  for  tests  on  sufficiently  soft  specimens  at  sufficiently  high  strain-rates, 
experimenters  should  abandon  the  practice  of  presenting  their  data  as  “stress-strain  curves”  and  simply 
report  the  measurements  that  are  valid.  For  the  torsional  Kolsky  bar  test  of  Nie  et  al.  [2013],  this 
consists  of  the  torque  history  at  the  load  cell  and  the  history  of  the  angle  of  twist  at  the  other  face  of  the 
specimen.  These  histories  should  be  provided  for  each  valid  test,  not  averaged  over  multiple  tests.  Such 
data,  in  conjunction  with  numerical  simulations  of  the  tests,  could  be  used  in  the  process  of  validating 
a  constitutive  model  that  has  already  been  calibrated.  In  this  regard,  optical  measurements  of  the  angle 
of  twist  at  the  end  of  the  incident  bar  adjacent  to  the  specimen  should  be  used  in  place  of  the  traditional 
strain  gage  measurements  on  the  incident  bar.  Unlike  the  analysis  of  strain  gage  data,  the  analysis  of 
optical  data  does  not  require  that  the  specimen  be  in  dynamic  equilibrium  or  that  the  incident  and  reflected 
pulses  do  not  overlap  at  the  measurement  location. 

The  above  procedure  could  not  be  used  to  actually  ccdibrate  the  parameters  in  a  viscoelastic  consti¬ 
tutive  model  unless  the  twist  and  torque  histories  are  part  of  a  larger  data  set  for  use  with  a  numerical 
inverse  method.  The  issue  here  is  nonuniqueness,  for  which  there  are  several  sources.  Even  if  a  numerical 
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inverse  method  were  used  to  optimize  the  fit  to  the  measured  stress  and  twist  histories,  it  seems  likely  that 
several  different  calibrations  might  fit  the  limited  data  equally  well,  particularly  in  view  of  the  variability 
that  exists  in  biological  specimens  and  the  fact  that  several  reasonable  criteria  for  optimality  could  be 
used.  Each  calibration  would  likely  result  in  a  different  numerical  prediction  of  the  strain  history  at  the 
load  cell.  Since  this  strain  history  cannot  be  measured  in  the  torsional  Kolsky  bar  test,  there  would  be  no 
criterion  for  selecting  one  of  these  calibrations  over  the  other.  Also,  significant  radial,  hoop,  and  axial 
stresses  exist  if  the  specimen  is  sufficiently  soft  and  the  strain-rate  sufficiently  high,  yet  these  normal 
stress  components  are  not  measured  in  the  torsional  Kolsky  bar  test.46  Furthermore,  the  measured  torque 
(or  the  torsional  shear  stress  Tqz )  may  be  insensitive  to  some  of  the  parameters  that  determine  these  normal 
stresses  as  was  demonstrated  for  the  case  of  a  Mooney-Rivlin  model.  Different  model  calibrations  might 
fit  the  twist  and  torque  histories  equally  well  or  nearly  as  well  and  yet  predict  drastically  different  normal 
stresses,  but  there  would  be  no  way  to  determine  which  (if  any)  of  these  predictions  are  correct. 

Appendix  A:  Details  of  numerical  methods 

One-dimensional  dynamic.  We  solve  the  wave  equation  (3-8)  along  with  boundary  conditions  (3-3) i 
and  (3-12)  numerically  using  a  finite  difference  scheme.  Initially,  an  explicit  second-order,  centered- 
time  and  centered-space  method  was  implemented.  This  method  produced  reasonable  results  but  had 
a  tendency  to  be  dispersive,  creating  erroneous  ripples  and  other  small  oscillatory  noise.  The  results 
presented  here  use  a  Crank-Nicolson  method  applied  to  the  centered-space  differences,  thus  making  it  a 
semi-implicit  method;  i.e.,  the  approximations  of  the  spatial  derivatives  were  taken  as  an  average  of  the 
values  for  the  current  time  step  and  the  next  time  step.  This  alteration  trades  the  dispersive  error  for  a 
slightly  dissipative  error.  The  length  of  the  specimen  (L  =  1.7  mm)  was  divided  up  into  equal  intervals 
(A z  =  0.0085  mm)  so  that  u(z,  t)  was  evaluated  at  200  points;  a  time  step  At  was  chosen  to  satisfy  the 
Courant  stability  criterion  for  the  explicit  method,  cAt/Az  <  1. 

Three-dimensional  quasistatic.  Meshes  were  generated  in  Cubit  (V12.1,  Sandia  National  Laboratory). 
Simulations  were  performed  using  Sierra  SolidMechanics  (Adagio  4.24,  Sandia  National  Laboratory). 
Adagio  is  an  implicit,  nonlinear  preconditioned  conjugate  gradient  solver.  Reduced  integration  on  HEX8 
element  meshes  was  performed.  Due  to  the  highly  constrained  nature  of  the  problem,  extreme  care  had 
to  be  taken  to  ensure  the  mesh  was  structured  and  symmetric  to  prevent  mesh  artifacts.  Postprocessing 
of  simulation  results  was  carried  out  in  ParaView  (V3.14.0,  Kitware)  and  MATLAB  (The  MathWorks). 

The  solutions  obtained  from  the  fully  three-dimensional  simulations  were  rotationally  symmetric  and 
enabled  a  reduction  of  the  data.  Variables  were  evaluated  along  five  radial  lines  at  a  single  angle  and 

1  1  T 

equally  spaced  axial  locations  (0  =  0  and  Z  =  0,  ^L,  ^L,  |L,  L).  This  postprocessing  was  performed  in 
ParaView  to  enable  further  analysis.  The  area  integrals  used  a  composite  Simpson’s  method  in  MATLAB. 
Similarly,  volumetric  averages  were  calculated  from  the  five  area  integrals  utilizing  a  five -point  Simpson 
rule.  Continuous  plots  over  the  radius  of  the  specimens  have  been  linearly  interpolated  between  elements. 

Three-dimensional  dynamic.  Several  demanding  simulation  requirements  arise  when  modeling  shear 
wave  propagation  in  soft,  nearly  incompressible  materials.  Figure  3  (for  the  one-dimensional  simulations) 

46  Measurements  of  axial  stress  in  quasistatic  torsion  tests  are  common,  and  such  measurements  could  perhaps  be  done  in 
torsional  Kolsky  bar  tests  as  well.  However,  it  seems  unlikely  that  radial  or  hoop  stresses  could  ever  be  measured. 
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shows  that,  when  a  specimen  is  far  from  dynamic  equilibrium,  large  gradients  in  strain  exist.  These  large 
gradients  cannot  be  captured  in  a  numerical  simulation  without  a  sufficiently  fine  mesh.  In  Simulation  /, 
there  were  changes  in  strain  of  0.4  over  distances  of  approximately  0.6  mm.  A  large  number  of  elements 
along  the  shortest  dimension  of  the  specimen  geometry  requires  extremely  large  meshes  to  model  the 
entire  specimen  in  a  three-dimensional  simulation.  A  large  reduction  in  the  number  of  elements  is  gained 
by  utilizing  the  rotational  symmetry  of  the  problem  and  effectively  changing  it  to  a  two-dimensional 
problem  by  simulating  a  thin  wedge  with  periodic  boundaries.  An  explicit  code  would  require  very  short 
time  steps  for  these  small  elements.  Simulations  were  performed  using  Adagio  with  HHT  integration 
[Hilber  et  al.  1977].  This  enabled  taking  larger  time  steps  while  maintaining  stability  and  introducing  a 
small  damping  to  higher  frequencies. 

The  software  used  for  meshing  and  postprocessing  was  the  same  as  for  the  quasistatic  simulations. 
Utilizing  the  rotational  symmetry  of  the  problem,  variables  were  evaluated  along  axial  lines  at  a  single 
angle  (0  =  0)  and  five  equally  spaced  radial  locations  from  Rj  to  R„.  For  the  simulations  with  a  vis¬ 
coelastic  specimen,  we  used  the  Viscoelastic  Swanson  Model  in  Sierra  SolidMechanics.  The  parameters 
in  the  elastic  paid  of  that  model  were  chosen  so  that  it  reduced  to  a  Mooney-Rivlin  model. 


Appendix  B:  Mathematical  formulation  of  a  smooth  loading  pulse 

The  simulations  of  the  torsional  Kolsky  bar  test  require  a  loading  pulse  (i.e.,  an  applied  average  strain- 
rate  history  ys)  that  ramps  up  smoothly  to  a  constant  value  before  subsequently  unloading  (Figure  2a). 
In  this  appendix,  we  formulate  a  parametrized  loading  pulse  using  a  Hermitian  smoothing  spline.  The 
time  integral  of  the  average  strain-rate  ys  can  then  be  used  to  obtain  an  applied  average  shear  strain  ys 
and  ultimately  the  applied  displacement  boundary  condition  (2-6). 

We  use  a  Hermitian  smoothing  spline  of  degree  5  as  the  basis  for  the  ramp-up  portion  of  the  loading 
pulse  (see  [Fitzpatrick  and  Scheidler  2013]): 

H5(x)  =  lO.r3  -  15x4  +  6.x5,  0<x<l.  (B-1) 


This  function  increases  from  0  at  x  =  0  to  I  at  x  =  I  and  has  zero  first  and  second  derivatives  at  both  of 
these  end  points.  Following  [Fitzpatrick  and  Scheidler  2013],  one  can  use  H${x)  to  construct  a  strain-rate 
function  that  ramps  up  from  zero  to  ymax  over  the  time  interval  [0,  tr  ],  maintains  this  constant  value  over 
the  time  interval  [tr,  tr  +  tp\,  and  then  unloads  to  zero  over  the  time  interval  \tr  +  tp,  td]: 


KsO)  = 


0, 

Ymax, 

Kmax[i-//5(f~(?;+^) 

0, 


t  <  0, 

0  <  t  <  tr, 

tr  <t  <  tr  +  tp , 
ri  T  tp<t  <  td , 
t  >  td. 


(B-2) 


Here  tr  is  referred  to  as  the  rise  time,  i.e.,  the  time  it  takes  to  reach  a  constant  strain-rate,  tp  is  the  duration 
of  the  constant  strain-rate  plateau,  and  td  is  the  total  duration  of  the  pulse.  In  the  simulations,  we  made 
the  additional  assumption  that  the  loading  pulse  is  symmetric;  that  is,  the  unloading  duration  td  —  (tr  +tp) 
equals  the  rise  time  so  that  td  =  2 tr  +  tp. 
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Note  that 


Ks(0)  =  Ys(tr)  —  Ys(tr  +  tp)  —  0,  Xs(0)  =  Ys(tr)  =  Ys(tr  +tp)  =  0.  (B-3) 


Since  ys  and  ys  are  identically  zero  during  the  constant  strain-rate  plateau  as  well  as  for  t  <  0  and  t  >  t,j, 
the  conditions  (B-3)  guarantee  that  ys  is  twice  continuously  differentiable.  This  degree  of  smoothness 
was  necessary  to  avoid  an  unwanted  ringing  in  the  simulations  that  would  otherwise  result  from  a  lower 
order  polynomial. 

An  explicit  expression  for  ys(t)  is  rather  messy,  and  in  practice,  it  is  simpler  to  work  with  the  integral 
of  H5  : 


H5(a) 


H5(y)dy  =  x4(x2 


3  A'  + 


(B-4) 


Integrating  each  case  in  (B-2)  with  respect  to  time,  assuming  a  symmetric  loading  pulse,  and  using 
H5(l)  =  \  gives 


ys(t)  = 


0, 

.  r 

Ymax 
Ymax 


\tr  +  (t  —  tr)\, 

\tr  +  (t  -  tr)  -  tr H5( 


t  ( tr~\~tp ) 


)] 


Kmax [ h  T  tp], 


t  <  0, 

0  <  t  <  tr, 
tr  <  t  <  tr  +  tp , 

tr  +  tp  <  t  <  2 tr  +  tp, 
t  >  1  tr  +  tp  =  td- 


(B-5) 


Equations  (B-2)  and  (B-5)  are  smooth  functions  that  qualitatively  reproduce  typical  strain  and  strain-rate 
loading  conditions  generated  in  Kolsky  bar  tests. 
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